NASA-CR— 199759 


V - 

The Telecommuhications and Data 
Acquisition Progress Report 42-123 

July-September 1995 

Joseph H. Yuen 

Editor 


(NASA— CR-199759) THE 
TELECOMMUNICATIONS ano Data 
ACQUISITION PRCGRFSS REPORT 92-127 
(JPL) 159 p 


•3 3/3, 


November 15, 1995 



National Aeronautics and 
Space Administration 

Jet Propulsion Laboratory 
California Institute of Technology 
Pasadena, California 


N 96- 16 664 
— THRU — 

N 9 6-16692 
Unc I as 


0 0 9 C 4 0 4 


The Telecommunications and Data 
Acquisition Progress Report 42-123 

July-September 1995 

Joseph H. Yuen 

Editor 


November 15, 1995 



National Aeronautics and 
Space Administration 

Jet Propulsion Laboratory 

California Institute of Technology 
Pasadena, California 


The research described in this publication was carried out by the Jet Propulsion Laboratory, California 
Institute of Technology, under a contract with the National Aeronautics and Space Administration. 

Reference herein to any specific commercial product, process, or service by trade name, trademark, 
manufacturer, or otherwise, does not constitute or imply its endorsement by the United States Government 
or the Jet Propulsion Laboratory, California Institute of Technology. 



Note From the Editor 


Since issue 42-121, published on May 15, 1995, The Telecommunications and 
Data Acquisition Progress Report has been available electronically to all its readers 
on the World Wide Web at http://tda.jpl.nasa.gov/progress_report. Printed copies 
were also produced but, as the Editor’s Note in that issue stated, the ultimate goal 
was to publish solely in electronic form. Now that goal is close to being realized. 
Beginning with issue 42-125, due to be published on May 15, 1996, The TDA 
Progress Report will be published entirely electronically at the above-mentioned 
URL. This issue and the next, 42-124, will be the last issues for which printed copies 
are produced; as a convenience for readers who are currently on our distribution 
list, however, we will distribute a hard copy of the table of contents for each issue. 
Readers with questions or concerns regarding this change are welcome to contact 
the editor. 



Preface 


This quarterly publication provides archival reports on developments in programs 
managed by JPL’s Telecommunications and Mission Operations Directorate (TMOD), 
which now includes the former Telecommunications and Data Acquisition (TDA) Office. 
In space communications, radio navigation, radio science, and ground-based radio and 
radar astronomy, it reports on activities of the Deep Space Network (DSN) in planning, 
supporting research and technology, implementation, and operations. Also included are 
standards activity at JPL for space data and information systems and reimbursable 
DSN work performed for other space agencies through NASA. The preceding work is 
all performed for NASA’s Office of Space Communications (OSC). 

TMOD also performs work funded by other NASA program offices through and 
with the cooperation of OSC. The first of these is the Orbital Debris Radar Program 
funded by the Office of Space Systems Development. It exists at Goldstone only and 
makes use of the planetary radar capability when the antennas are configured as science 
instruments making direct observations of the planets, their satellites, and asteroids of 
our solar system. The Office of Space Sciences funds the data reduction and science 
analyses of data obtained by the Goldstone Solar System Radar. The antennas at all 
three complexes are also configured for radio astronomy research and, as such, conduct 
experiments funded by the National Science Foundation in the U.S. and other agencies 
at the overseas complexes. These experiments are either in microwave spectroscopy or 
very long baseline interferometry. 

Finally, tasks funded under the JPL Director’s Discretionary Fund and the Caltech 
President’s Fund that involve TMOD are included. 

This and each succeeding issue of The Telecommunications and Data Acquisition 
Progress Report will present material in some, but not necessarily all, of the aforemen- 
tioned programs. 
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Sensitivity of Planetary Cruise Navigation 
to Earth Orientation Calibration Errors 

J. A. Estefan 

Navigation and Flight Mechanics Section 
W. M. Folkner 

Tracking Systems and Applications Section 


A detailed analysis was conducted to determine the sensitivity of spacecraft navi- 
gation errors to the accuracy and timeliness of Earth orientation calibrations . Anal- 
yses based on simulated X-band (8.4-GHz) Doppler and ranging measurements ac- 
quired during the interplanetary cruise segment of the Mars Pathfinder heliocentric 
trajectory were completed for the nominal trajectory design and for an alternative 
trajectory with a longer transit time. Several error models were developed to char- 
acterize the effect of Earth orientation on navigational accuracy based on current 
and anticipated Deep Space Network calibration strategies. The navigational sensi- 
tivity of Mars Pathfinder to calibration errors in Earth orientation was computed for 
each candidate calibration strategy with the Earth orientation parameters included 
as estimated parameters in the navigation solution. In these cases, the calibration 
errors contributed 23 to 58 percent of the total navigation error budget, depend- 
ing on the calibration strategy being assessed. Navigation sensitivity calculations 
were also performed for cases in which Earth orientation calibration errors were not 
adjusted in the navigation solution. In these cases , Earth orientation calibration 
errors contributed from 26 to as much as 227 percent of the total navigation error 
budget. The final analysis suggests that . not only is the method used to calibrate 
Earth orientation vitally important for precision navigation of Mars Pathfinder, but 
perhaps equally important is the method for inclusion of the calibration errors in 
the navigation solutions. 


I. Introduction 

Radio metric data, particularly two-way coherent Doppler and range, have been used to navigate 
robotic spacecraft since the inception of planetary exploration. For a spacecraft in interplanetary cruise 
or transit, much of the information content inherent in the data for position determination comes from 
the signature imposed on the station-spacecraft radio signal by the Earth’s rotation [1-3]. The diurnal 
signature in the radio metric data yields information about the right ascension and declination of the 
spacecraft with respect to the direction of the Earth’s spin axis at the time of observation. The orientation 
of the Earth, as a function of time, must be known with respect to inertial space in order to effectively 
utilize the radio metric data to deduce spacecraft position with respect to the target planet. Errors in 
Earth orientation thus lead to targeting errors for spacecraft approaching other planetary bodies. 
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Evidence of the need to adequately account for Earth orientation errors came as early as April 1965 
when flight project navigation teams for the Rangers VII and VIII lunar probes observed a large difference 
in station longitude solutions for all deep-space stations using radio metric data [4]. This was later 
determined to be the result of improper Earth orientation calibration. As a result, Earth orientation 
calibration methods were later refined to support the Mariners IV and V planetary exploration missions. 

To assess the effect of Earth orientation calibration errors on interplanetary cruise navigation for both 
current and future Deep Space Network (DSN) Earth orientation calibration techniques, a navigation error 
analysis of the Mars Pathfinder approach scenario was performed. Mars Pathfinder has the most stringent 
planetary cruise navigation requirements of any currently planned mission. Other Mars lander missions 
similar to Pathfinder are being studied. Navigation performance for these future missions may exhibit 
different sensitivity characteristics to Earth orientation calibration errors since the sensitivity is trajectory 
dependent. In this study, two Mars approach trajectories were evaluated, the nominal Pathfinder cruise 
trajectory with arrival at Mars on July 4, 1997, and a second trajectory with a longer transit time. Clearly, 
restricting the study to only two trajectories is far from encompassing the entire range of possible planetary 
approach scenarios. Moreover, actual navigation performance will vary depending on targeting point and 
targeting requirements, the data type and arc length, filtering strategy, and observation geometry, which 
could vary depending on each launch opportunity (especially spacecraft right ascension and declination at 
encounter). The study of two “representative” trajectories, while limited, at least provides some insight 
into the possible range of navigational uncertainties caused by Earth orientation calibration errors. 

Other types of navigation problems have varying sensitivity to Earth orientation error. For spacecraft 
in close orbit about another planetary body, such as Magellan or Mars Global Surveyor (MGS), the 
primary signature on the spacecraft radio signal is imposed by the orbit about the planet. Doppler 
measurements can be used to determine all spacecraft orbital elements in most cases, and the resultant 
orbit determination is largely insensitive to Earth orientation errors. If, however, the orbit determination 
using Doppler data is not accurate enough to meet the mission requirements, two-station differenced- 
Doppler or narrow-band very long baseline interferometry (VLBI) observations can be used for improved 
orbit determination in some cases. In fact, the Magellan project utilized differenced-Doppler data for this 
purpose. A detailed sensitivity analysis of differenced-Doppler navigation to Earth orientation calibration 
errors is not presented in this article, but a cursory approximation is given in Appendix A. In contrast to 
spacecraft in close planetary orbit, the Galileo and Cassini spacecraft will be in long-period (~1 20-day) 
orbits about Jupiter and Saturn, respectively. These represent an intermediate case between planetary 
approach and low planetary orbit, so some sensitivity to Earth orientation errors might be expected. 
Onboard optical images of planetary satellites will be an important data type in determining the orbits 
for Galileo and Cassini. The added complexity of blending onboard optical data with radio metric data 
precluded this study from assessing the navigation sensitivity to Earth orientation errors for these outer 
planet orbiters. 

In this article, a navigation error analysis is described that was used to assess the impact of various 
Earth orientation calibration strategics on predicted spacecraft orbit determination accuracies during 
interplanetary cruise. Section II provides the fundamental framework for defining the principal param- 
eters that arc used to characterize Earth orientation, while Section III focuses on the Earth orientation 
calibration process used by the DSN. These discussions are followed by a description in Section IV of 
the origin and format of the functional requirements levied on the DSN tracking system by the flight 
projects. In Section V, a simple information content analysis is presented to obtain a rough estimate of 
the influence of Earth orientation errors on Doppler cruise navigation performance. Section VI describes 
the assumptions used in a linear covariance analysis to evaluate the sensitivity of spacecraft navigational 
accuracies to Earth orientation calibration errors for two Mars Pathfinder approach scenarios. Various 
Earth orientation calibration strategies art; described, together with tracking data simulation and error 
modeling assumptions. Results and key observations from the numerical assessment are summarized and 
discussed at the conclusion of the article. 
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II. Earth Orientation Parameters 

The Earth is an oblate, spinning body that undergoes precession and nutation due to the torques 
exerted upon it by the Sun, Moon, and other planets. The north pole of a body-fixed (crust-fixed) 
coordinate system varies unpredictably with respect to the spin direction, due to internal dynamics of the 
Earth and its atmosphere (a process called “polar motion”). Similar effects cause the Earth’s rotation 
rate to vary unpredictably. (The variations in the rotation rate are several times larger than the polar 
motion variations.) 

The orientation of a body in inertial space can be completely described by three Euler angles. Because 
the Earth rotates rapidly, the three angles describing the orientation of the surface with respect to inertial 
space vary rapidly with time. Conventionally, the orientation of the Earth is described by five angles that 
vary slowly with time, rather than by three rapidly varying angles. These five angles are described in 
greater detail below. 

In the development of the 1980 International Astronomical Union (IAU) theory of nutation [5], the 
concept of the celestial ephemeris pole (CEP) was introduced. The CEP was defined such that there 
are no nearly diurnal motions of the CEP with respect to either space-fixed (inertial) or body-fixed, 
coordinates. For a rigid body with no polar motion, the CEP corresponds to the body axis about which 
the body is spinning. 

The motion of the CEP in space-fixed coordinates, due to precession and nutation, can be described by 
the two angles, 0 and £, where e is the obliquity (inclination of the equatorial plane to the Earth’s orbital 
plane) and 0 is the intersection of the equator and orbit with respect to a fixed equinox. The variation in 
the Earth's rotation about the CEP affects the time at which celestial objects cross the apparent meridian 
and is measured by a quantity called Universal Time (UT) (specifically, Universal Time 1, or “UT1”). 
Variations of the CEP in body-fixed coordinates are measured by the quantities polar motion “X” and 
polar motion “Y”. 

Because of the random variation of UT1 and polar motion (along with imperfect modeling of precession 
and nutation), an accurate description of Earth orientation requires continual monitoring. VLBI data 
can be used to determine all components of Earth orientation with 5 nr ad or better accuracy (1-sigma). 1 
Because VLBI measurements require correlation of large volumes of data from ground stations separated 
by large distances, there is usually a time delay between the acquisition of raw VLBI data and the 
processing of the data that determines the Earth orientation angles. This processing delay is currently 2 
to 3 days for DSN VLBI measurements made for rapid determination of Earth orientation (i.e., TEMPO 
measurements, described in Section III); the delay is longer for VLBI data from external services. Satellite 
laser ranging (SLR) or Global Positioning System (GPS) data can also be used to determine polar motion 
and small changes in UT1 with shorter data-pro cessing times but are not able to directly measure all 
five Earth orientation angles. Atmospheric angular momentum (AAM) data are highly correlated with 
variation in UT1 and the length of the day (LOD), a parameter proportional to the rate of change of 
UT1. Therefore, AAM data, both measurements and forecasts, have been used to improve predictions 
for both TJT1 and LOD [6]. 

Precession/nutation models with parameters adjusted to fit the observed space-fixed motion of the 
CEP, e.g. ,[7,8] , have an accuracy of 5 nrad or better over the time of the fit. These models can be used 
to predict precession and nutation for periods of about 1 year before discrepancies systematically exceed 
5 nrad. Figure 1 shows a comparison of the daily correction calibrations (from 1991 to 1995) with a 
model by Steppe et al. fit to data through the end of 1993 [9]. The nutation corrections are with respect 


1 Earth orientation accuracies are often quoted in a variety of units. An angle of 5 nrad is approximately equal to 1 milliarc- 
second (mas). An angular rotation of 5 nrad corresponds to a change in position on the surface of the Earth (equatorial 
displacement) of about 3 cm. A change in UT1 of 1 millisecond (ms) corresponds to an angle of about 15 mas, which is 
equivalent to an angle of 75 nrad, or to an equatorial displacement of roughly 50 cm. 
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Fig. 1. A comparison of a nutation correction model to observations. The correction angles 
(a) sin(e) 8y and (b) 8e are corrections to the IAU (1976) nutation model [5]. (A change of 8e or 
sin(e) 5u i of 1 mas corresponds to a shift in the inertial position of a point on the Earth's surface 
of about 3 cm.) 


to the 1976 IAU precession model [10,11] and 1980 IAU nutation model [5]. The corrections are currently 
about 10 mas (~50 nrad) and are increasing with time. It can be seen that the predictions of the model 
fit to data through the end of 1993 agree with the later measurements to an accuracy of about 1 mas for 
about 1 year. 

UT1 and polar motion (collectively referred to as UTPM throughout this article) vary randomly due 
to the interaction of the atmosphere and the crust. UT1 varies much more rapidly than polar motion. 
Random variation in UT1 can be characterized by an integrated random walk, while polar motion behaves 
approximately as an integrated Gauss-Markov process [12]. UT1 varies by an amount corresponding to 
an angle of 1 mas in about 1 day, so continual, rapid calibration is required to be able to completely 
describe Earth orientation to 1-mas accuracy. 


III. Earth Orientation Calibrations 

At, the Jet Propulsion Laboratory (JPL), Earth orientation calibrations are currently determined by 
the DSN’s Time and Earth Motion Precision Observations (TEMPO) activity. TEMPO, which became 
operational in late 1983, was chartered to provide an operational Earth orientation service both to support 
JPL’s spacecraft navigation efforts and to serve the worldwide community [13]. TEMPO supports Earth 
orientation calibration by performing VLBI measurements at regular intervals (currently twice per week) 
using the DSN’s 70-m antenna subnetwork. (Prior to 1983, Earth orientation calibrations were provided 
by the DSN’s Tracking System Analytic Calibration (TSAC) activity, which produced calibrations based 
on monthly estimates of UTPM disseminated by the Bureau Internationale de l’Heure (BIH) in Paris, 
France.) 2 The Kalman Earth Orientation Filter (KEOF) is used to combine the TEMPO measurements 


2 T. F. Range, personal communication, Tracking Systems and Applications Section, Jet Propulsion Laboratory, Pasadena, 
California, February 1995. 
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with other sources of Earth orientation information. By performing regular VLBI measurements and 
including AAM measurements and forecasts, together with other data from Earth orientation services 
outside JPL, the DSN can deliver, at any time, an Earth orientation calibration accurate to 50 nrad 
(1-sigma) [6]. Earth orientation accuracy is better for times 15 days or more in the past, for which a 
greater amount of processed VLBI data from external services is available. Earth orientation predictions 
are also delivered, with accuracies that degrade with time due to the random behavior of UTPM. Efforts 
are ongoing to improve these predictions both by modeling improvements within the KEOF and by better 
utilization of geodetic and AAM data. 3 

The standard DSN Earth orientation calibration file (referred to as a UTPM STOIC file) is a text file of 
polynomial coefficients that provides UTPM calibrations for 37 specified times. 4 Precession and nutation 
calibrations are not included. The limitation to 37 calibration times has implications for the accuracy of 
Earth orientation available to the end-user (e.g., navigation teams) because of the integrated-random walk 
characteristic of UTPM. Several flight projects utilize calibration files that span a year or more, giving 
10-day spacing (or more) between calibration times. Midway between respective calibrations at 10-day 
intervals, the expected (1-sigma) error in UT1 is about 0.4 ms (~20 cm) even if the calibration is perfect at 
the calibration times. This limitation, together with the lack of precession/nutation calibrations, has led 
to a new DSN calibration file— the Earth-Orientation Parameter (EOP) file — which includes precession 
and nutation corrections and has no limit on the number of calibration times. 5 6 7 It should be noted that 
all timing calibrations and their rates are given with respect to a reference time defined by atomic clocks, 
specifically, International Atomic Time (TAI). 


IV. Functional Requirements 

Navigation-related requirements for current and future missions are defined primarily by flight projects 
and future mission study teams. These requirements serve as a starting point to establish DSN ground 
support requirements to satisfy mission navigation. In the past, navigation requirements for calibrations 
such as Earth orientation, station locations, and transmission media typically have been arrived at in 
an ad hoc manner without thorough analysis. This practice has at times resulted in confusion and later 
cancellation of implementation plans to develop calibrations for which there was an erroneously believed 
need. 

Arguably, flight projects and future mission study teams find it more economical to simply adopt past 
calibration performance or to adopt anticipated improvements in the calibrations rather than conduct a 
parametric study in which all possible navigation calibrations are investigated. In order to meet mission 
navigation needs, the DSN has documented UTPM calibration capabilities and requirements for the 
tracking and navigation subsystems, 6,7 The current UTPM requirements are stated as: “(a) 30 cm (1- 


3 J. O. Dickey, personal communication, Tracking Systems and Applications Section, Jet Propulsion Laboratory, Pasadena, 
California, August 1995. 

4 In the early 1970s, all UTPM calibration data for mission operations were supplied in a single computer card deck called 
a PLATO deck (Platform Observables). The PLATO system replaced the former Timing and Polynomial (TPOLY) 
computer program for generating separate timing calibration data [14], For contingency purposes, a smaller and simpler 
backup program was developed to generate PLATO-style decks that could be delivered rapidly in the event PLATO was 
not operable. This program was called STOIC (Standby Timing Operation In Contingencies) — hence, the frequently 
encountered convention “STOIC” file or, more appropriately, “UTPM STOIC” file. Sometimes, these files are referred to 
by their historical convention as “TPOLY” files or simply as “TP” arrays. 

7 DSN Tracking System Interfaces , Earth Orientation Parameter Data Interface (TRK-2-21), DSN System Requirements 
Detailed Interface Design , JPL 820-13, Rev. A (internal document), Jet Propulsion Laboratory, Pasadena, California, 
April 19, 1985. 

6 DSN System Functional Requirements and Design: Tracking System (1988 Through 1993), JPL 821-19, Rev. C (internal 
document), Jet Propulsion Laboratory, Pasadena, California, pp. 3-20, April 15, 1993. 

7 NOCC Subsystem Functional Requirements: Navigation Subsystem (1988 Through 1993), JPL 822-18, Rev. A (internal 
document), Jet Propulsion Laboratory, Pasadena, California, pp. 3-7-3-8, May 15, 1988. 
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sigma) in each component, predictive, for the days on which the calibrations are generated; (b) 5 cm 
(1-sigma) in each component, non-predictive, for periods through 14 days prior to the day on which the 
calibrations are generated; (c) 5 to 25 cm (1-sigma) in each component of polar motion, non-predictive, 
for periods from 1962 through 1984; and (d) 10 to 40 cm (1-sigma) in UT1, non-predictive, for periods 
from 1962 through 1984. ” 8 

The exact origin of the 30-cm real-time knowledge requirement is not widely known, although it 
is clear that it was arrived at via the common practice of synthesizing past flight project navigation 
team requirements and what the current calibration activity claimed could be delivered in terms of 
accuracy and timeliness. There is a common misconception that the 30-cm functional requirement for 
all three components of UTPM was driven by Magellan mission requirements. In actuality, the Magellan 
30-cm requirement was inherited directly from the Galileo project for a 30-cm real-time UTPM knowledge 
requirement. 9 The UTPM requirements levied by future missions (e.g., Cassini, MGS) vary from flight 
project to flight project and are subject to change. Therefore, mission-specific requirements will not be 
presented here. It is fair to state that an effort is under way to update the overall Earth orientation 
calibration functional requirements for Mars Pathfinder (precession/nutation as well as UTPM) based on 
the analysis presented in this article. 


V. Information Content Analysis 

Early analytic studies suggest that Earth orientation uncertainties result in equivalent uncertainties in 
the instantaneous location of tracking stations, which leads to a degradation in the apparent quality of the 
radio metric data used for navigation [15-17]. As noted in the introductory remarks, timing (UT1) errors 
in particular can lead to an erroneous prediction of the spacecraft coordinates near planetary encounter. 
Much of Doppler data’s information content, when acquired during interplanetary cruise, comes from 
the diurnal signature of the Earth’s rotation. This is evident in a simple analytic representation of the 
instantaneous range rate, p, observed by an Earth-based tracking station [1—3] : 

p = v r 4- r s uj e cos<5sina; e £ + ( — Aa + A A 4- uj e AUTl)r s uj e cos S cos uj e t (1) 


Here, v r denotes the spacecraft radial velocity with respect to the Earth; r s is the distance of the station 
from the Earth’s spin axis, u> e denotes the rotation rate of the Earth, and t is measured from the nominal 
time the spacecraft crosses the tracking station’s meridian. The 6 is the instantaneous declination of the 
spacecraft, Aa the correction to the a priori value of spacecraft right ascension, A A the correction to the 
station longitude, and AUT1 the correction to rotation about the spin axis. There are, of course, other 
parameters to be estimated. Moreover, this simple model neglects the additional geometric strength that 
comes from the motion of the Earth about the Sun and the use of multiple tracking stations. Nevertheless, 
this model is useful to illustrate (to first order) the effect of Earth orientation errors on the Doppler data. 

It is clear from Eq. (1) that an error in rotation about the Earth’s spin axis would directly affect the 
right ascension estimate. For example, at Deep Space Station (DSS) latitudes (~35 deg), a 1-ms timing 
error is equivalent to a longitude error of about 0.4 m, or a right ascension error of about 0.07 prad 
[13]. Polar motion affects the spacecraft position estimate by producing displacements in the station 
spin radius, longitude, and height above the equator. These displacements can be as large as 10 m 
if not properly calibrated. Equation (1) expresses the spacecraft right ascension and declination with 
respect to the Earth’s equator of date. Errors in precession and nutation models can lead to errors in 
the transformation of the “of-datc” right ascension and declination estimate into the inertial coordinate 

8 Ibid. 

9 S. N. Molian and W. L. Sjogren, “Revised Navigation Requirement Specification for the VRM Mission Requirements Docu- 
ment 630-6 and Preliminary Spacecraft Instrumentation Requirements Document (SIRD),” JPL Interoffice Memorandum 
314. 10-348, Rev. 1 (internal document), Jet Propulsion Laboratory, Pasadena, California, September 22, 1983. 
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system of the planetary ephemeris. Precession/nutation modeling errors are rarely significant for Earth- 
orbiting spacecraft, where the observational data tie the spacecraft orbit much more tightly to a local 
coordinate system. For interplanetary spacecraft, the trajectory determined by Earth-based radio metric 
data must be related to the position of a distant planet, 

VI. Navigation Error Analysis 

To investigate the effect of various levels of Earth orientation calibration accuracy on interplanetary 
cruise navigation, an error covariance analysis of the Mars Pathfinder approach segment was performed. 
The Mars Pathfinder approach scenario was selected because it has the most stringent planetary ap- 
proach navigation requirements of any currently planned mission. Future Mars lander missions may 
utilize different trajectory designs and potentially could exhibit a lesser or greater level of sensitivity to 
Earth orientation calibration errors than those presented herein. This analysis is intended to serve as a 
representative model. 

A. Calibration Strategies (Test Cases) 

In order to study the effect of Earth orientation calibrations on Mars Pathfinder cruise navigation, 
six test cases were developed to cover a wide range of possible Earth orientation calibration strategies. 
The level of calibration errors, which are a function of time, depends on the amount of data included in 
creation of the calibration files and on the timeliness of their deliveries. Precession/nutation calibrations 
were not included in this study since it is possible to predict the corrections for about a year with an 
accuracy approaching ~1 mas. All cases of Earth orientation studied here assume a basic set of VLBI 
measurements that can provide this level of precession/nutation accuracy (cf., Section II). 

The Earth orientation calibration cases are characterized by the uncertainty in UT1 and polar motion 
as a function of time and by the correlations between the errors at different times. The reference day for 
the Earth orientation calibration cases is the day on which the navigation solution is performed. Figure 2 
shows the assumed uncertainty in UT1 for the six cases, while Fig. 3 illustrates the assumed polar motion 
uncertainties. To simplify the analysis, the indicated level of polar motion uncertainty was assumed for 
both the X and Y components independently, even though actual measurements show that uncertainty 
in predictions for polar motion Y increase about 30 percent slower than for polar motion X [12]. Table 1 
gives additional information about the statistics of each Earth orientation case. In all cases, the errors 
due to the potentially sparse array of calibration times imposed by the STOIC file have been neglected 
since this effect can be removed either by use of the EOP file or by use of a STOIC file that spans the 
shortest time possible. 

The baseline Earth orientation case is the current nominal DSN capability and is indicated as TEMPO 
in the figures and in Table 1 [6]. 10 This is based on two DSN VLBI measurements per week combined 
with data available from other sources. It is assumed that the last processed TEMPO VLBI measure- 
ment of UT1 was acquired 5 days before the KEOF filter run and that the KEOF filter run is performed 
1 day prior to the navigation solution. For UT1 prediction, the rate of change in UT1 is important. 
The UT1 rate is dependent on the last two processed TEMPO measurements. For this particular case, 
UT1 was characterized as a first-order Gauss-Markov random process with a 1-sigma steady-state un- 
certainty of 0.11 ms and a 5-day correlation time until 7 days prior to the navigation solution; from this 
time forward, UT1 uncertainty was characterized by an integrated random walk (through the time of the 
navigation solution). The current KEOF filter solutions include the TEMPO VLBI measurements and 
AAM measurements and forecasts to give the stated capability for UT1 accuracy on any given day. In 
addition, daily VLBI measurements from external services are included in the KEOF filter solutions to 
provide the steady-state uncertainty of 0.11 ms for times in the past. (The 0.11 ms is larger than the 


10 A. P. Freedman, “Polar Motion Prediction With KEOF,” JPL Interoffice Memorandum 335.2-92.01 (internal document), 
Jet Propulsion Laboratory, Pasadena, California, March 5, 1992. 
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Fig. 2. 

Modeled UT1 errors versus time for six different calibration strategies. 
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Fig. 3. Modeled polar motion errors versus time for three different calibration 
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Table 1. Earth orientation calibration accuracies 
for various strategies. 


TAI - UT1 Polar motion 


Calibration strategy 

I -sigma, 

ms 

1-sigma, 

mas 

-21 days 

0 days 

—21 days 

0 days 

TEMPO (current,) 

0.11 

0.71 

1.6 

5.6 

Delayed TEMPO 

0.11 

1.32 

1.6 

7.7 

Delayed TEMPO - AAM 

0.11 

1.93 

1.6 

7.7 

TEMPO + CPS (CODE) 

0. 11 

0.23 

1.6 

2.3 

CPS (CODE) + 2 VLBI/mo 

0.18 

0.39 

1.6 

2.3 

GPS (JPL) -f 2 VLBI/mo 

0.18 

0.25 

1.6 

2.3 
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quoted measurement uncertainties in order to accommodate possible offsets and drifts between various 
Earth orientation services.) In addition, the KEOF includes polar motion determinations from SLR 
measurements of Earth-orbiting satellites, which are available up to 5 days before the filter run (6 days 
before the navigation solution). The uncertainty in each component of polar motion was modeled as 
first-order Markov with a 5-day correlation time and a 1-sigma steady-state uncertainty of 1.6 mas, up 
to 6 days before the navigation solution, at which time the uncertainty was modeled as a random walk 
increasing to the final time. For times later than 6 days before the navigation solution, polar motion was 
modeled as a random walk that approximated the observed polar motion statistics [12]. (An integrated 
Gauss-Markov process was not used due to the difficulty in implementing it in the covariance analysis 
software.) 

To investigate the importance of timely Earth orientation calibration delivery, two cases were included 
with a 6-day delay between the KEOF filter run and the navigation solution. The case labeled “de- 
layed TEMPO” in Figs. 2 and 3 and in Table 1 is identical to the baseline case except for an additional 
6-day delay. For the delayed TEMPO case, UT1 uncertainty is assumed to grow as an integrated ran- 
dom walk 13 days before the navigation solution, and polar motion begins to grow as a random walk 
11 days before the navigation solution. A third case, labeled “delayed TEMPO - AAM,” is identical 
to the delayed TEMPO case except that AAM data in the KEOF solution are not included. Without 
the AAM data, the UT1 uncertainty begins growing as an integrated random walk 13 days prior to the 
navigation solution, but at a faster rate. The polar motion uncertainty is identical for those two cases, 
i.e., delayed TEMPO and delayed TEMPO - AAM. 

Measurements of GPS satellites have been used extensively for geodetic purposes, and GPS data have 
a demonstrated capability to measure polar motion and LOD. (Recall LOD is directly related to UT1 
rate.) For the past 2 years, the Center for Orbit Determination, Europe (CODE) has been producing 
daily measurements of polar motion and LOD with a week or more delay between data acquisition and 
Earth orientation delivery. There are plans for the DSN to begin rapid processing of GPS data to sup- 
plement and partially replace TEMPO VLBI measurements in an effort to reduce loading on the DSN’s 
70-m subnetwork. If the measurements of LOD are uncorrelated (i.e., “white”), then including LOD 
measurements implies a random walk noise on UT1. Three cases of VLBI and GPS data combinations 
are included here assuming that the GPS LOD measurements are uncorrelated. The actual noise char- 
acteristics are under investigation. If the LOD measurements turn out to be correlated, 11 then the effect 
of GPS Earth orientation calibrations on navigation error may be different than the results presented in 
this article. 

The fourth Earth orientation case, labeled “TEMPO + GPS(CODE),” assumes the current level of 
external VLBI measurements, the current two TEMPO VLBI passes each week, plus GPS polar motion 
and LOD measurements with a 1-day processing time. For this case, UT1 was assumed to behave as a 
Gauss-Markov process with a 5-day correlation time and a 1-sigma steady-state uncertainty of 0.11 ms 
until 7 days before the navigation solution, at which time the UT1 uncertainty was assumed to grow as 
a random walk at a level characteristic of the CODE GPS LOD deliveries. Each component of polar 
motion is described by a Gauss-Markov process with a 5-aay con elation time and a 1.6- mas steady- 
state uncertainty (1-sigma) until 2 days before the navigation solution, at which time the polar motion 
uncertainty increases as a random walk. (This polar motion uncertainty model is assumed for all three 
cases that include GPS data.) 

Because the current DSN plan is to utilize GPS LOD measurements so as to acquire fewer VLBI 
measurements, and because the number of external VLBI services has been steadily declining, the fifth 
Earth orientation case assumes that only two VLBI measurements are acquired per month and combined 


11 Preliminary studies of JPL GPS-derived LOD measurements exhibit “nonwhite” behavior on time scales longer than 3-5 
days, A. P. Freedman, personal communication, Tracking Systems and Applications Section, Jet Propulsion Laboratory, 
Pasadena, California, August 18, 1995. 
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with GPS measurements. The case, labeled “GPS (CODE) + 2 VLBI /mo” in the figures and in Table 1 
assumes a 10-day delay in processing the VLBI measurements. This delay may occur as a result of 
reducing the load on the 70-m subnetwork, whereby VLBI measurements are acquired using the 34-m 
subnetwork. This strategy would require tapes to be used to record the VLBI data and shipped back to 
JPL for processing. With a 10-day delay between VLBI data acquisition and final processing, the latest 
VLBI measurement the KEOF could possibly include would be 11 days before the navigation solution 
with a worst-case delivery of 25 days before the navigation solution. For this case, it was assumed that 
the UT1 uncertainty behaved as a Gauss-Markov process with a 5-day correlation time and a 1-sigma 
steady-state uncertainty of 0.18 ms until 18 days prior to the navigation solution. This higher steady- 
state uncertainty is due to the lack of daily VLBI measurements from external services and reflects the 
uncertainty from using daily GPS LOD measurements to interpolate between VLBI UT1 measurements. 
At 18 days before the navigation solution, the UT1 uncertainty is assumed to grow as a random walk at 
a level characteristic of the CODE LOD measurements. 

The current DSN plan is to have in place a JPL rapid GPS processing system for Earth orientation. 
The 3-year implementation cycle will begin in fiscal year 1996 with provisional operations beginning 
as early as fiscal year 1998. This will give way to a fully operational system by fiscal year 1999. 12 
The processing implementation plan is under development but, as a test, there have been daily GPS 
solutions for LOD performed since late 1994. These solutions do not span a long enough time period 
to provide a good statistical measure of performance, but preliminary results indicate the JPL LOD 
measurements may be twice as accurate as the CODE deliveries. The sixth Earth orientation case, 
labeled “GPS(JPL) + 2VLBI/mo,” is identical to the previous test case, GPS(CODE) + 2VLBI/mo, 
except that at 18 days before the navigation solution, the uncertainty in UT1 is assumed to begin a 
random walk behavior with a slower growth rate. 

B. Mars Pathfinder Tracking and Error Modeling Assumptions 

The Mars Pathfinder spacecraft will directly enter the Martian atmosphere from Earth transfer orbit 
for landing on the Martian surface. Other missions (e.g., Cassini, MGS) will either fly by the target 
planet or enter orbit through a series of orbital correction maneuvers. The primary atmospheric entry 
constraint for Mars Pathfinder is the flight path angle, the angle between the incoming velocity vector 
of the spacecraft and the vector normal to the Martian atmosphere. If this angle is too large (shallow), 
the spacecraft may overheat before parachute deployment, and if the angle is too small (steep), excess 
pressure may develop that could potentially damage the spacecraft’s aeroshell from ablation. This entry 
angle constraint is expected to place the most stringent requirements on calibration of Earth orientation. 
A secondary requirement is to target the spacecraft to land within a predetermined landing footprint on 
the Martian surface. The size of the landing footprint is 100 km x 300 km. 

The Mars Pathfinder spacecraft will be spin stabilized throughout its interplanetary cruise to Mars 
and will communicate through its high-gain antenna. The onboard telecommunications system has an 
X-band (7.2-GHz) uplink/X-band (8.4-GHz) downlink radio system, which will be used to acquire Doppler 
and ranging measurements and to transmit science and engineering telemetry data. The nominal launch 
window is a 30-day launch period beginning on December 5, 1996. 13 Arrival at Mars is scheduled to 
occur on July 4, 1997. The launch vehicle will be targeted so that it will not impact the Martian sur- 
face. In the first 60 days after launch, two trajectory correction maneuvers (TCMs) will be performed 
to remove the effects of launch vehicle injection errors and to remove the targeting bias. A third TCM 
(TCM-3) is scheduled to be executed 60 days prior to Mars atmospheric entry. The critical navigation 
event time is just before the final maneuver (TCM-4). Five days prior to TCM-4, a navigation solu- 
tion will be generated to design the final maneuver. The maneuver design command parameters will 

12 S. M. Lichten, personal communication, Tracking Systems and Applications Section, Jet Propulsion Laboratory, Pasadena, 
California, January 1905. 

M At the time this article went to print, the actual launch window was not yet fixed since the mission profile and spacecraft 
launch mass were still being refined. 
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be uplinked to the spacecraft for execution from 10 to 15 days before atmospheric entry. Expected tra- 
jectory uncertainties for this critical navigation delivery have been carefully studied by Thurman and 
Kallemeyn 14, 15,16 via linear covariance analysis and Monte Carlo simulation. The covariance analysis 
assumptions adopted herein to assess the sensitivity of the critical Mars Pathfinder navigation solution 
to various Earth orientation calibration strategies were derived in large part from these earlier navigation 
performance assessments. 

The nominal Mars Pathfinder trajectory is a so-called “Type I” trajectory, where the heliocentric 
longitude of the spacecraft changes by less than 180 deg between launch and arrival. An alternative 
“Type II” trajectory, where the heliocentric longitude of the spacecraft changes by more than 180 deg 
and less than 360 deg between launch and arrival, was originally considered for Mars Pathfinder. Analysts 
who first studied the Type II trajectory option suggest that the principal reasons the Type I trajectory 
option was preferred were (1) to attempt to minimize 70-m antenna conflicts between Mars Pathfinder at 
arrival and the Galileo mission at Jupiter, (2) to shorten the cruise time from ~11 months to ~7 months, 
which would yield less consumables in terms of propellant, and (3) to attain a more favorable geometry 
for the spacecraft to remain at Earth-point during cruise while maximizing the Sun’s exposure to the 
solar arrays. 17 (The Sun-probe-Earth angle is small for this mission.) A navigation error analysis for the 
Type II option was included in this assessment because some future missions to Mars (including MGS) 
will utilize Type II trajectories. 


Table 2. Assumed data arc lengths for Mars Pathfinder 
navigation analysis. 


Trajectory 

Launch/arrival date 

Data arc specification 

Begin, days 3 

End, days b 

Length, days 

Type I 

January 3, 1997/ 

L + 60 

M- 15 

107 


July 4, 1997 (fixed) 




Type II 

December 2, 1996/ 

L + 236 

M— 15 

107 


November 10, 1997 (fixed) 





a L = launch. 
b M = Mars arrival. 


The tracking data arcs assumed for the covariance analysis are shown in Table 2 for both the Type I 
and Type II trajectories. X-band two-way coherent Doppler and ranging data were simulated over these 
intervals. DSN coverage varied according to the nominal DSN data acquisition schedule specified in the 
Mars Pathfinder Navigation Plan. 18 For the Type I transfer phase (L-b60 days to M — 45 days), the DSN 
coverage wets taken to be one 4-h pass/week per complex; during the Mars approach phase (M — 45 days to 
Mars arrival), continuous coverage was assumed; and for the TCM-3 phase, one 8-h pass/day (continuous 
for 12 h before and after TCM-3) was assumed over the interval TOM ± 3 days. The same data arc 


14 S. W. Thurman, “Orbit Determination Filter and Modeling Assumptions for MESUR Pathfinder Guidance and Navigation 
Analysis,” JPL Interoffice Memorandum 314.3-1075 (internal document), Jet Propulsion Laboratory, Pasadena, California, 
October 15, 1993. 

15 Navigation Plan: Preliminary Version , Pathfinder Flight Project, JPL D-11349 (internal document), Jet Propulsion 
Laboratory, Pasadena, California, December 1993. 

16 Navigation Plan: Critical Design Review Version , Mars Pathfinder Project, JPL D-11349 (internal document), Jet Propul- 
sion Laboratory, Pasadena, California, July 1994. 

17 V. M. Pollmeier, personal communication, Navigation and Flight Mechanics Section, Jet Propulsion Laboratory, Pasadena, 
California, March 1995. 

18 Navigation Plan: Critical Design Review Version , op. cit. 
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length was used for the Type II trajectory; thus, simulated data points began at L 4- 236 days. In an 
effort to minimize the effects of potential station or complex outages while maximizing the angle-finding 
capability of the ranging data, tracking passes were scheduled to alternate between DSN complexes. 

The Doppler and ranging data were assumed to have measurement uncertainties of 0.09 mm/s (60 s 
average) and 2 m, respectively (1-sigma). Although recent X-band Doppler data residuals are typically 
smaller than 0.09 mm/s, a higher Doppler uncertainty was assumed in order to reflect the low-frequency 
power of the solar plasma noise spectrum that is not properly characterized by the root-mean-square of 
the residuals. 19 A 20-min integration time was assumed for each data point (for both data types). 

The Mars Pathfinder trajectories were integrated from initial position and velocity conditions (epoch 
state) using models for the dynamic forces on the spacecraft. The modeled gravitational forces were due 
to the masses of the Sun and the planets; relative locations of these bodies were based on the JPL DE200 
ephemeris. Other forces modeled were nongravitational accelerations due to solar radiation pressure 
(SRP), gas leaks from valves and pressurized tanks, and attitude maintenance activity. In addition, 
TCM-3 maneuver execution errors were modeled. 

Parameters estimated by the data reduction algorithm (a variant of the sequential Kalman filter [18]) 
included a wide array of dynamic and observational error sources categorized as (1) spacecraft epoch 
state, (2) spacecraft nongravitational force modeling errors, (3) maneuver execution errors, (4) errors 
in the orbital elements of the Earth and Mars, (5) systematic Doppler and ranging error biases, (6) 
transmission-media zenith delay calibration errors for the ionosphere and troposphere, (7) crust-fixed 
station location errors, and (8) Earth orientation calibration errors for UTPM. All of these error sources 
and their assumed a priori and steady-state values are summarized in Table 3. A priori uncertainties 
for the spacecraft initial state were large enough to leave it essentially unconstrained, while nongravita- 
tional forces were modeled as first-order Gauss-Markov random processes. (Note that all nongravitational 
forces except the slowly varying SRP accelerations were modeled using a stochastic gas leak model and 
are lumped under the category a NGA” in the table, where NGA denotes nongravitational accelerations.) 
TCM-3 execution errors (for the TCM-4 delivery) were modeled as random biases in all three body-fixed 
components. The uncertainty in the Earth-Mars ephemeris was taken from the JPL DE234 ephemeris 
error covariance by Standish, 20 but constrained with the knowledge that the orientation of the Earth s 
orbit is now known to 15 nrad [19]. For processing the two-way ranging data, the filter model included a 
bias parameter associated with each ranging pass from each station in order to approximate the slowly 
varying nongeometric delays in the ranging measurements that are caused principally by station delay 
calibration errors and uncalibrated solar plasma effects. The spacecraft spin rate, detectable in the 
Doppler signature, was estimated as a Gauss-Markov process with a 5-day correlation time. Uncertainty 
in knowledge of the station locations was assumed to be 10 cm for each component. This station location 
uncertainty is expected to be characteristic of the new DSN beam-waveguide (BWG) antennas once sur- 
veys are complete. More accurate station locations exist for antennas for which VLBI data are available, 
including the 70-m antennas and the 34-m high-efficiency (HEF) antennas. 21 

Although this study is restricted to the orbit determination problem and does not address the influence 
of guidance errors on navigational accuracy, it is important to note that, upon completion of TCM-4, 
the contribution of maneuver execution errors to the overall guidance dispersions are expected to be 
negligible. This was demonstrated in preflight error analyses and is discussed in greater detail in the 
Mars Pathfinder Navigation Plan. 22 

,5) W. M. Folkner, “Effect of Uucalibrated Charged Particles on Doppler Tracking,” JPL Interoffice Memoi andum 

335.1-1)1-005 {internal document), Jet Propulsion Laboratory, Pasadena, California, March 1, 1994. 

20 K. M. Standish, “The JPL Planetary Ephemerides, DE234/LE234,” JPL Interoffice Memorandum 314.6-1348 (internal 

document), Jet Propulsion Laboratory, Pasadena, California, October 8, 1991. 

21 In actuality, the Mars Pathfinder spacecraft will be “uplink-limited” and will, therefore, require use of the 34-m II Eh 

antennas for telecommunication. A more conservative assessment is made herein by assuming the 34-m BWG antennas. 

22 Navigation Plan: Critical Design Review Version , op. cit. 
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Table 3. A priori and steady-state uncertainties for orbit 
determination error model parameters. 


Estimated parameter set 

Uncertainty, 1 a 

Remarks 

Spacecraft epoch state 

A priori 

Constant parameters 

Position components 

100 km 


Velocity components 

1 m/s 


Nongravitational force model 

Solar radiation pressure (SRP) 

Steady-state 

First-order Markov 

Radial ( G r ) 

5% of nominal 

60-day correlation time 

Transverse ( G x /G y ) 

5% of nominal 


Gas leaks (NGA) 

Steady-state 

First-order Markov 

Radial (a r ) 

2 x 10" 12 km/s 2 

5-day correlation time 

Transverse ( a x /a y ) 

2 x 10" 12 km/s 2 

5-day correlation time 

Maneuver execution error model 

TCM-3 (AV*,AV yi AV Z ) 

A priori 

Constant parameters 

(for TCM-4 delivery) 

10~ 2 m/s 


Planetary ephemerides error model 

Earth- Mars ephemeris 

A priori 

Constant parameters 

Orbit orientation (3 Euler angles) 

15 nr ad 


Longitude with respect to periapsis 

10 nrad 


Semimajor axis (A a/a) 

5 parts in 10 1 1 


Eccentricity (Ae) 

3 parts in 10 10 


Ground system error model 

Range biases 

A priori 

Constant parameters 

(one per station per pass) 

1 m 


Transponder bias 

Steady-state 

First-order Markov 

(ranging data only) 

1 m 

0.5-day correlation time 

Doppler spin bias 

Steady-state 

First-order Markov 

(Doppler data only) 

10 -2 mm/s 

5-day correlation time 

Transmission media 

Steady-state 

First-order Markov 

Zenith troposphere 

5 cm 

0.1-day correlation time 

Zenith ionosphere 

5 x 10 16 e/m 2 

0.2-day correlation time 

DSN station coordinates 

A priori 

Constant parameters 

(crust-fixed r s , z/ l5 A) 

10 cm 

(uncorrelated) 

Earth orientation 

(cf. , Section VI. A) 

(cf., Section VI. A) 

Timing (UT1) 

Polar motion (X,Y) 


C. Encounter Geometry 

Because much of the strength of the Doppler and ranging data comes from the signature imposed by 
the rotation of the Earth, interpretation of the covariance analysis results is aided by understanding the 
encounter geometry. 

The spacecraft position in the Earth-spacecraft direction is directly measured by ranging data. Doppler 
data help determine the other two components of the spacecraft position, which lie in the plane of the sky. 
There is a well-known weakness in determining spacecraft declination from Doppler data for spacecraft 
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near zero declination [1-3]. Spacecraft declination can be inferred from ranging data using tracking 
stations located in both northern and southern latitudes [20]. Figure 4 shows the Pathfinder Type I 
trajectory on the plane of the sky as viewed from Earth. As seen in the figure, encounter occurs near 
zero declination. Because of this encounter geometry, the spacecraft declination will probably depend 
upon ranging data, and the declination uncertainty should exhibit sensitivity to station delay calibration 
errors. In contrast, the Type II trajectory has a relatively large, negative encounter declination, as shown 
in Fig. 5. For the Type II encounter, the Doppler data will have a larger role in determining spacecraft 
declination, which should thus be less sensitive to station delay calibration errors. 
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Fig. 4. Mars Pathfinder Type I trajectory as viewed from Earth; 
shown at 5-day intervals from the beginning of the data arc to 
encounter. 
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Fig. 5. Mars Pathfinder Type II trajectory as viewed from Earth; 
shown at 5-day intervals from the beginning of the data arc to 
encounter. 


A typical uncertainty ellipsoid for the spacecraft position on approach would have principal axes ap- 
proximately aligned with the plane-of-sky axes, with a much smaller uncertainty in the Earth-Mars 
direction than in the other two components (assuming ranging data are included). Planetary approach 
trajectories are typically described in aiming plane (B-plane) coordinates. 23 Figure 6(a) shows the rela- 
tionship of the D - plane components to the plane-of-sky components for the Pathfinder TypeJ encounter. 
The approach direction, S, is nearly parallel to the radial (Earth-Mars) direction, r. The -R direction is 
in the plane normal to the approach direction, S, and approximately parallel to the direction of increasing 
declination, 8, while the -T direction is in the plane normal to the approach direction and approximately 


23 For a complete description of the B-plane coordinate system, please refer to Appendix B. 
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(a) 


-T 


Fig. 6. The S-plane components for Mars Pathfinder approach tra- 
jectories with respect to plane-of-sky coordinates: (a) Type I and 
(b) Type II. 

parallel to the direction of decreasing right ascension, -a. For the Type I trajectory, the well-determined 
component is approximately in the direction of approach. A small position uncertainty in this direction 
is expressed in the 5-plane system as a small uncertainty in the time from encounter, i.e., linearized 
time of flight (LTOF). The position uncertainty is approximately related to the LTOF uncertainty by the 
approach velocity. For the Type I trajectory, the approach velocity is about 5.5 km/s (a 1-s uncertainty 
in LTOF corresponds to a position error of about 5.5 km). An error in right ascension, such as might be 
caused by a UT1 calibration error, will appear in the B • T component. 

Figure 6(b) shows the relationship of the 5-plane components to the plane-of-sky components for the 
Type II trajectory. The direction of the spacecraft approach to Mars, — S, is about 11 deg from the 
direction of decreasing right ascension, -5. The -R direction is about 23 deg from the declination axis, 
5. The — T direction is about 23 deg from the Earth-Mars direction, r. For this trajectory, an error in 
right ascension will be reflected mainly in LTOF. For the Type II trajectory, the approach velocity is 
approximately 3.9 km/s. 

Mars Pathfinder navigation is required to deliver, prior to the final maneuver (TCM-4), a trajectory 
estimate with less than a 1-percent probability of exceeding the entry angle requirement. The latest 
assessment of the Type I flight path entry angle requirement is ±1 deg (99 percent), which implies a 
requirement on the navigation delivery corresponding to a 3-sigma uncertainty of 21 km in the magnitude 
of the impact parameter. 24 Stated another way, the entry corridor is 42-km wide, as depicted in Fig. 7. 

D. Resuits 

In the covariance studies performed, a careful model was constructed for the time-dependent Earth 
orientation errors shown in Figs. 2 and 3. This model would be somewhat difficult to implement into 
the operational Orbit Determination Program (ODP), 25 which currently does not have a statistical reset 
capability or an integrated random walk model such as the one used in this analysis. Because of this 
limitation, the effect of each Earth orientation calibration strategy on the total orbit determination error 
was calculated in two ways. For the first estimation method, the contribution to orbit determination 
error from Earth orientation was determined with UT1 and polar motion included (i.e., estimated) in 
the navigation solution, with correctly modeled time-dependent a priori uncertainties. In the second 
estimation method, the contribution to orbit determination error was assessed under the assumption that 
the Earth orientation calibration errors were ignored (i.e., not estimated) in the navigation solution. 


24 P. K. Kallemeyn, personal communication, Navigation and Flight Mechanics Section, Jet Propulsion Laboratory, 
Pasadena, California, March 1995. 

25 The ODP is a large institutional software system used for research and navigation support of flight operations, N. D. Pana- 
giotacopulos, J. W. Zielenbach, and R. W. Duesing, An Introduction to JPLs Orbit Determination Program , JPL 1846-37 
(internal document), Jet Propulsion Laboratory, Pasadena, California, May 21, 1974. 
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Fig. 7. Mars Pathfinder entry corridor and landing accuracy 
requirements (99 percent). (Note: the Earth equatorial plane is nearly 
parallel to the direction to Mars.) 


Figure 8 shows the contribution to the total navigational uncertainty for the nominal (Type I) Mars 
Pathfinder trajectory from all error sources described in Section VI.B except Earth orientation. The 
covariance analysis results given below are expressed in £?-plane components referred to the Earth mean 
equator of J2000 (EME2000), as described in Appendix B. Contributions were computed in a manner such 
that the sum of each error source, when added in quadrature, gives the total navigation uncertainty in a 
root-sum-square sense [21]. The critical navigation parameter for Pathfinder approach is the magnitude 
of the impact parameter, denoted |B|. Recall from the previous discussion that |B| is related to the 
flight path entry angle. For the nominal Pathfinder Type I approach trajectory, the B • T uncertainty, 
denoted C B is nearly the same as the uncertainty in |B|. In general, the relationship of the component 
uncertainties cr B g,a B and C| B | depends upon the choice of the targeted entry point. 

The uncertainties in arrival time (LTOF) are very small because of the approach direction nearly 
coinciding with the Earth-Mars direction, which is well determined by ranging data. The scale for 
LTOF in Fig. 8(c) is 3 s and corresponds to a position uncertainty of about 15 km. The major error 
source (other than Earth orientation) for B • T (~right ascension) is the anomalous nongravitational 
accelerations (NGAs). The B • R (^declination) uncertainty has roughly equal contributions from data 
noise, nongravitational forces, and station delay calibrations for ranging data. The 1-m accuracy of the 
range bias calibrations assumed for the covariance analysis has been inferred from observations of the day- 
to-day consistency of Mars Observer ranging data residuals [22]. This assumption should be interpreted 
cautiously since the systematic effects in the Mars Observer range biases could have been absorbed by 
other spacecraft trajectory parameters, such as nongravitational accelerations. Fortunately, this is not 
an issue for Mars Pathfinder since the critical navigation component, |B|, is almost entirely in the right 
ascension direction. Further, this navigation error analysis was not intended to be the “official” Mars 
Pathfinder analysis. The principal purpose here was to provide a quantitative measure of the relative 
importance of potential error sources, specifically, Earth orientation calibration errors. 

Figure 9 illustrates the contribution to the total orbit determination uncertainty from each case of 
Earth orientation calibration error described in Section VI. A. Here, it is seen that Earth orientation 
calibration errors are a significant source of error for Mars Pathfinder in the critical B ■ T and |B| 
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Fig. 8. Relative contributions of the principal error sources (other than Earth orientation) to the total Mars 
Pathfinder orbit determination uncertainty^ Uncertainties are shown in B-plape coordinates with respect to 
the mean Earth equator of 2000: (a) B • R (-declination) uncertainty, (b) B • T (-right ascension) uncertainty, 
(c) LTOF (time of encounter) uncertainty, and (d) uncertainty in the magnitude of the impact parameter, IBI. 


components . 26 Earth orientation calibration error is less significant in the B • R and LTOF components. 
The lack of sensitivity to Earth orientation in the LTOF direction is due to the fact that the approach 
direction is nearly aligned with the Earth-Mars direction; therefore, LTOF is well determined by the 
ranging data. The spacecraft declination (nearly aligned with the B ♦ R direction) is determined largely 
by ranging data at northern and southern latitude stations since, at the low encounter declination, the 
Doppler data do not contribute much to the determination of declination. Because the declination is 


26 Recall that errors due to precession and nutation were neglected from this analysis; thus, the formal Earth orientation 
calibration errors are strictly due to UTPM calibration errors. 
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Fig. 9. Relative contributions of Earth orientation to the total Mars Pathfinder orbit determination uncertainty. 
Uncertainties are shown in 8-plane coordinates with respect to the mean Earth equator of 2000: (a) B • R 
(-declination) uncertainty, (b) B • T (-right ascension) uncertainty, (c) LTOF (time of encounter) uncertainty, and (d) 
uncertainty in the magnitude of the impact parameter, IBI. Uncertainties are given for both the case where Earth 
orientation parameters were estimated in the navigation solutions and for cases where Earth orientation was not 
adjusted in the navigation solution. 


determined principally by ranging data with an assumed accuracy of 1 m, there is not much sensitivity 
to Earth orientation errors for the calibration strategies studied here, which all give Earth orientation 
errors smaller than 1 m at the Earth’s surface. 

In the case of current DSN Earth orientation calibration performance, assuming a delivery of the cali- 
bration files on the day of the critical navigation solution from TEMPO VLBI data, Fig. 9(d) shows that 
Earth orientation errors contribute approximately 39 percent of the 1-sigma |B| (flight path entry angle) 
requirement of 7 km for the case where UTPM parameters were included in the navigation solution, 
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and 64 percent of the allowable error if UTPM were ignored in the navigation solution. The two Earth 
orientation calibration cases with delayed delivery show contributions to navigation uncertainty that are 
significantly larger. The delayed calibration cases are most likely unacceptable for the Mars Pathfinder 
mission. The optimistic Earth orientation case, in which the current twice-weekly TEMPO VLBI mea- 
surements are augmented with daily GPS data, shows a much smaller contribution to the navigation 
uncertainty than the nominal TEMPO case. The two calibration strategies with daily GPS data com- 
bined with reduced VLBI observations (2 VLBI/month) are comparable to the nominal TEMPO case. 
In contrast to the nominal TEMPO case, the GPS-based calibrations exhibit smaller differences between 
the strategy of including UTPM parameters and statistics in the navigation solution and the strategy of 
ignoring the UTPM parameters in the navigation solution. 

Figure 9 shows a reduced sensitivity to Earth orientation errors when the UTPM parameters are 
estimated, along with the trajectory parameters, in the navigation filter. This improvement is large for 
the cases with poorest UTPM accuracy. The improvement is coupled to the assumptions about the level 
of nongravitational forces affecting the spacecraft. If there were no nongravitational forces acting on the 
spacecraft, or if the nongravitational forces were perfectly known, then the spacecraft would provide a 
reference against which Earth orientation changes could be measured using Doppler data. If there were 
large nongravitational forces affecting the spacecraft that are not well known, then the spacecraft could not 
be used as a reference against which Earth orientation changes could be measured. Because the Pathfinder 
spacecraft will be a simple, spinning platform, the nongravitational forces affecting it are assumed here 
to be well modeled. Because of this assumption, when the Earth orientation uncertainties increase 
beyond a certain level, the navigation filter begins to rely on the assumed level of nongravitational force 
uncertainties and can improve upon the a priori knowledge assumed for Earth orientation parameters. 
This would not be true for a spacecraft with larger uncertainties in the nongravitational force model. 

Because of the different encounter geometry, the covariance analysis results for the Type II trajectory 
are quite different from the Type I trajectory. No attempt was made to quantify the critical navigation 
component, |B|, since the Type II trajectory will not be used for Mars Pathfinder and the choice of the 
targeted point for this study was arbitrary. The B-plane component ^uncertainties should be interpreted 
in such a manner that the critical component could be more like B • R or B T, depending on the choice 
of the targeted point. 

Figure 10 shows the navigation uncertainty from all error sources for the Type II trajectory with the 
exception of Earth orientation calibration error. The LTOF uncertainty is about a factor of six larger for 
the Type II case than for the Type I case because the approach direction is not aligned with the Earth- 
Mars direction. The scale in Fig. 10(c), 3 s, corresponds to a position uncertainty of about 12 km due to 
the approach velocity of 3.9 km/s. Nongravitational forces, Mars ephemeris uncertainty, and data noise 
are seen to be the dominant sources of error (other than Earth orientation) for the other components. 
The B • T component is most closely aligned with the Earth-Mars direction at encounter and, hence, is 
the best determined component. The uncertainty in B R (^declination) shown in Fig. 10(a) for the 
Type II trajectory is less sensitive to ranging calibration errors and more sensitive to station location 
errors than is the Type I case. This is a reflection of the large, negative encounter declination enabling 
Doppler data to influence the determination of declination. 

Figure 11 shows the contribution of Earth jprientation calibration errors to the orbit determination un- 
certainty for the Type II trajectory. The B • R uncertainty is much more dependent on Earth orientation 
than is the Type I case. This sensitivity is related to the determination of declination by the Doppler 
data, which are sensitive to Earth platform errors. The sensitivity of declination to Earth orientation can 
be seen to be principally due to polar motion errors since cases with identical polar motion uncertainties, 
but different UT1 uncertainties, have the same effect on the B • R uncertainty. The B ■ T compo- 
nent shows some sensitivity to UT1 errors since the T direction is mostly in the Earth-Mars direction 
but partly in the direction of increasing right ascension. UT1 errors have a larger effect on LTOF since the 
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Fig. 10. Relative contributions ot the principal error sources (other than Earth orientation) to the total 
orbit determination uncertainty for a Mars Pathfinder Type II approach scenario.JJncertainties are shown in 
B-plan^coordinates with respect to the mean Earth equator of 2000: (a) B • R (-declination) uncertainty, 
(b) B • T (-right ascension) uncertainty, (c) LTOF (time of encounter) uncertainty, and (d) uncertainty in the 
magnitude of the impact parameter, IBI. 


approach direction is more closely aligned with the right ascension direction. The nominal TEMPO Earth 
orientation errors would be one of the larger sources for error in B • T and a moderate source of error in 
B • R for this trajectory. The GPS-based cases contribute less to the navigation uncertainty in B • R and 
B • T than the nominal TEMPO case, but result in large errors in LTOF. 
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Fig. 11. Relative contributions of Earth orientation for a Mars Type II approach trajectory. Uncertainties are 
shown in B-plane coordinates with respect to the mean Earth equator of 2000: (a) B • R (-declination) 

uncertainty, (b) B • T (-right ascension) uncertainty, (c) LTOF (time of encounter) uncertainty, and (d) uncertainty 
in the magnitude of the impact parameter. IBI. Uncertainties are given for both the case where Earth orientation 
parameters were estimated in the navigation solutions and for cases where Earth orientation was not adjusted in 
the navigation solution. 


VII. Summary and Conclusions 

A numerical assessment measuring the sensitivity of spacecraft delivery errors to the accuracy and 
timeliness of Earth orientation calibrations was completed for two interplanetary cruise scenarios derived 
from the Mars Pathfinder mission set. This study was motivated by the fact that, to date, errors in Earth 
orientation (i.e., precession/nutation, polar motion, and variation in Earth rotation rate) are still capable 
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of contributing significantly to the composition of the noise signature on radio metric data acquired by 
the DSN. These errors can thus lead to degraded spacecraft navigational accuracies if not adequately 
calibrated. 

Results from the navigation sensitivity analysis concurred with the expected outcome that not only 
is Earth orientation calibration performance important in determining spacecraft navigational accuracy, 
but so is the timeliness of the calibration file deliveries. Based on the analyses presented in this article, 
the current best DSN Earth orientation calibration performance provided by the TEMPO activity yielded 
a contribution of about 39 to 64 percent of the total navigation error budget for the critical component 
of the nominal Mars Pathfinder Type I trajectory, depending on the navigation filtering strategy being 
used. These results assumed line-of-sight data types (i.e., two-way Doppler and range) were used in 
the navigation process. Use of differential data types could reduce the sensitivity to Earth orientation 
calibration errors. 

Variations on the current DSN calibration method representing delayed TEMPO deliveries of the 
calibration files as well as delayed deliveries without use of AAM data were also assessed. Results for 
these cases showed that Earth orientation calibration errors dominated the total navigation error budget, 
irrespective of the trajectory type. Furthermore, a very large penalty was paid when the Earth orientation 
parameters were not adjusted in the navigation solution. 

With the advent of GPS-based ground observations as a viable Earth orientation calibration system and 
the ongoing effort to reduce the loading on the DSN 70-m subnetwork, new Earth orientation calibration 
techniques are being devised. Statistical models representing examples of these calibration strategies 
were constructed and their effect on the Mars Pathfinder navigation delivery error assessed. In the 
(optimistic) case where the current level of TEMPO calibrations (2 per week) was used in concert with 
daily GPS-based calibrations, the influence of UTPM calibration errors on overall navigation performance 
was, as expected, minimal. Under the current environment where there is continual pressure to reduce 
the number of DSN-based VLBI observations (again, addressing the 70-m antenna loading issue), this 
calibration strategy will probably not be attainable operationally. 

A sensitivity analysis was also performed for an operationally more realistic Earth orientation calibra- 
tion strategy in which GPS-based calibrations were used as the principal means of generating frequent 
(daily) Earth orientation calibration information, augmented with periodic VLBI-based measurements 
(~2 per month). (The GPS system alone cannot determine all components of Earth orientation and, 
thus, requires an external calibration source such as VLBI.) In this assessment, analysis results suggest 
that the contribution of UTPM errors to the total navigation error budget for the critical component 
of the nominal Mars Pathfinder trajectory lies somewhere between 43 and 55 percent, depending on the 
accuracy of the GPS deliveries. These results assumed that the UTPM parameters were adjusted in the 
navigation solution. The true level of accuracy will depend, of course, on the actual system implemented. 

Since the GPS calibration system is in the early stages of development, the statistical characteristics 
of the calibrations are not yet well determined. With the noise levels assumed for this analysis, the 
GPS-based Earth orientation calibrations appear to offer an advantage over the current TEMPO-based 
calibrations in that they relax the need for the navigation process to properly model the time-varying 
behavior of UTPM calibration errors. In addition, the proposed system is designed to provide rapid 
processing and timely deliveries of the calibration files to the flight projects. The overall performance 
(accuracy) levels, as evidenced in this study, were at or near the same level as the current DSN capability, 
perhaps only slightly better in some cases. 
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Appendix A 

Sensitivity of Planetary Orbiter Navigation to 
Earth Orientation — A Case Study for 
Differential Data Types 


For a spacecraft in orbit about another planet, Doppler data can be used to determine all components 
of its orbit except for a few particular geometries [23]. The accuracy with which the orbit is determined 
by means of Earth-based Doppler tracking depends upon several factors, including data accuracy and the 
accuracy of the spacecraft force models, particularly those due to the planet’s gravitational field. Using 
the Magellan radar mapping mission of Venus as an example, the uncertainty of the gravity field was 
such that the expected orbit uncertainty during the prime mission (for daily orbit solutions) was about 
15 km with two-way Doppler tracking alone [24]. Mars Global Surveyor (MGS) plans to achieve much 
better accuracy by solving for an improved gravity field based on an initial data set. The MGS strategy 
could not have been utilized for Magellan since the gravity field of Venus could not be sampled with a 
few weeks of radio metric data because of Venus’ slow rotation rate. 

The orbit determination accuracy achievable with Earth-based Doppler tracking in a two-way coherent 
mode is very insensitive to Earth orientation errors since the dominant signature in the Doppler data is 
due to the orbit of the spacecraft about the planet. This is in contrast to planetary approach navigation, 
where much of the information content in Doppler tracking data is influenced by the Earth’s rotation. 
To first order, Doppler data are insensitive to a rotation about the line of sight from the Earth-based 
tracking station to the spacecraft. The rotation about the line of sight can be determined by changes in 
the geometry due to the relative orbits of the Earth and the target planet about which the spacecraft is 
orbiting. Rotation about the line of sight is measured by the node angle, Q, with respect to the plane of 
the sky, which is defined in Fig. A-l as the plane normal to the Earth-spacecraft direction. In the figure, 
the orbit inclination, 2 , is the angle between the normal to the spacecraft orbit and the Earth -spacecraft 
direction; the line of nodes is the intersection of the orbit plane and the plane of the sky; the node with 
respect to the plane of the sky, Q, is measured in the plane of the sky from a reference direction to the line 
of nodes; and the argument of periapsis, oo, is the angle between the line of nodes and periapsis measured 
in the orbit plane. 


r (TOWARD EARTH) 



Fig. A-1 . Planetary orbiter geometry. 
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At times, the desired orbit accuracy is greater than what can achieved with two-way Doppler tracking 
alone. In the case of Magellan, the desired orbit accuracy was about 1 km for purposes of aligning 
radar images. This level of accuracy was better than what could be achieved using Doppler data alone. 
Differential data types such as differenced one-way Doppler (DOD), delta-differenced one-way Doppler 
(ADOD), or two-way minus three-way Doppler (2DM3D) have been shown to improve orbit determination 
accuracy [26]; the latter was used successfully for Magellan operations [27]. These differential data 
types are sensitive to Earth orientation errors. An assessment of the characteristic sensitivity to Earth 
orientation errors for planetary orbiter navigation when using differential data types is described below. 
The planetary orbiter scenario is based on the radar mapping phase of the Magellan mission. 

Consider the geometry drawn schematically in Fig. A-2. The plane of the figure is taken to be the 
Earth’s equatorial plane. (Note that Venus need not lie in the equatorial plane for the analysis to be 
valid.) Two stations at different complexes are located at the ends of the baseline vector with equatorial 
projection, b e . For illustration purposes, an orbit is considered with the orbit plane perpendicular to 
the Earth’s equatorial plane and with the normal to the orbit plane perpendicular to the Earth- Venus 
direction. 



Fig. A-2. Differenced Doppler measurement geometry used in the case study. 


DOD measurements are formed by differencing the one-w T ay Doppler signals received by two tracking 
stations separated by large distances [26]. These measurements give the difference in spacecraft line-of- 
sight velocity as observed by the two stations. (The 2DM3D measurements exhibit the same information 
content except for a slight difference resulting from use of a DSN uplink signal rather than the spacecraft 
onboard oscillator as the reference frequency.) For spacecraft at interplanetary distances, the DOD 
observable can be approximated as 


DOD “ s ( b ■ ;) “ r [ b ' ( v - f 7) + (; * “•) ■ b ] 


(A-l) 


where b is the baseline vector between the two tracking stations, r is the vector from the center of the 
Earth to the spacecraft with magnitude r, r is the rate of change of distance between the Earth and 
the spacecraft, v is the spacecraft velocity vector with magnitude u, and co e is the Earth’s rotation rate 
vector. By considering a DOD measurement for this special case, at the time when the spacecraft velocity 
is parallel to the Earth’s pole, the DOD observable can be further approximated as 


DOD « -v y b e cos H H — v z b z + u e b e cos H 
r r 


(A-2) 


where v y is the component of the spacecraft velocity in the equatorial plane and perpendicular to the 
Earth- Venus direction (and nominally zero at the measurement time), b e is the equatorial baseline length, 
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v z is the component of the spacecraft velocity parallel to the Earth’s pole, b z is the length of the projection 
of the baseline length onto the pole direction, and H is the hour angle between the baseline and the 
spacecraft. A rotation of the orbit about the Earth- Venus line by an angle 6D changes v y from zero to 
vSQ, which is directly observable in the DOD measurement. If the measurement occurred earlier (or later) 
in the orbit, where the spacecraft velocity vector was along the spacecraft-Earth direction, there would 
be no change in spacecraft velocity for a change in the orbit node. In this case, the DOD measurement 
would not be useful. The importance of performing differenced Doppler measurements at optimum times 
has been well documented in the literature (see, e.g., [25]). 

An error in UT1 introduces a bias in the hour angle, H , and, hence, in the DOD measurement. 
This can affect the determination of the spacecraft node angle. A change in the measurement due to a 
calibration error, 6UT 1, is approximately given by 


6 DOD « - b e sin H6 UT 1 


(A-3) 


This change will cause an error to be inferred in the rotation about the line of sight by an amount 


Sfl « -ujI tan H 6UT1 (A-4) 

V 

For DSN baselines (Goldstone-Madrid and Goldstone-Canberra), H can vary from about —30 to +30 
deg, outside of which the spacecraft will fall below the horizon of one of the complexes. DSN baselines 
have a mean equatorial length of about 8000 km. For the worst case where H = 30 deg, a 1-ms error 
in UT1 will bias the DOD measurements by about 0.02 mm/s. For an orbiter characteristic of Magellan 
during its mapping phase, with an average orbital velocity, v of about 5.5 km/s, and a line-of-sight 
distance of 1 AU, a 1-ms timing error in UT1 would lead to a node error of up to 0.08 mrad. With a 
semimajor axis of 10,000 km, this corresponds to an orbit error of about 0.8 km. (Since this is comparable 
to the desired orbit accuracy for Magellan, it was necessary to have UT1 calibrated with submillisecond 
accuracy in order to support the generation of daily orbit determination solutions.) 

In general, the maximum sensitivity of the differenced Doppler data to Earth orientation errors is of 
nearly the same magnitude as for the special case studied here. Sensitivity to Earth orientation errors 
can be an order of magnitude smaller if the data are acquired at times where the baseline hour angle is 
near zero and the spacecraft velocity at that time is in a direction where the data are sensitive to the 
spacecraft node. The size of the orbit errors also depends on (among a number of other factors) the shape 
of the orbit, the uncertainty in the gravity field, and the amount of Doppler data to be used in the “fit” 
(i.e., the data filtering process). 
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Appendix B 

Definition of Aiming Plane (B-Plane) Coordinates 


Planetary approach trajectories are typically described in aiming plane coordinates, often referred to 
as “5-plane” coordinates (see Fig. B-l). This coordinate system was originally conceived to simplify the 
targeting of a hyperbolic flyby trajectory and is defined by three orthogonal unit vectors, S, T, and_R, 
with the system origin taken to be the gravitational center of mass of the target planet [27]. The S is 
directed parallel to the incoming spacecraft asymptotic velocity vector relative to the target planet, while 
T is normally specified to lie in either the ecliptic plane (the mean plane of the Earth’s orbit) or the 
equatorial plane of the target planet J 27 In addition, T is directed perpendicular to S. The unit vector R 
completes an orthogonal triad with S and T, thus, R = S x T. 

The aim point for a planetary encounter is defined by the impact parameter, B, which approximates 
where the point of closest approach would be if the target planet had no mass and did notjieflect the flight 
path. The impact parameter B is directed perpendicular to S; therefore, it lies in the T - R plane. To 
gain insight into targeting accuracy, orbit determination errors are often characterized by the 1-sigma or 
3-sigma uncertainty in the respective “miss components” of B, namely, B • R and B T. These quantities 
are analogous to elevation and azimuth when specifying the impact point for terrestrial targets. 

The time from encounter is defined by the linearized time of flight (LTOF), a quantity which is a 
measure of the “time-to-go” from the current spacecraft position to the intersection of its asymptotic 
flight path and the aiming plane. LTOF provides a convenient time-to-go parameter because LTOF is 
not affected by changes in the B • R and B ■ T miss components. 28 Orbit determination errors are also 
characterized by the 1-sigma or 3-sigma uncertainty in LTOF. 

In lieu of using B • R and B • T uncertainties to measure targeting accuracy, a 1-sigma or 3-sigma 
5-plane dispersion ellipse (also shown in Fig. B-l) is often used. The semimajor (SMAA) and semiminor 
(SMIA) axes of the dispersion ellipse are related in quadrature to the uncertainties of B R and B • T. 
The angle Of gives the angle clockwise from T to the SMAA. 


TARGET PLANET 



Fig. B-1. The aiming plane (B-plane) coordinate system. 


27 For the analysis presented in this article, T was specified to lie in the Earth’s equatorial plane. 

28 R. A. Jacobson, “Linearized-Time-of-Flight Revisited,” JPL Engineering Memorandum 391-680 (Revised) (internal doc- 
ument), Jet Propulsion Laboratory, Pasadena, California, September 22, 1975. 
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Wind data measured during a field experiment were used to verify the analytical 
model of wind gusts . Good coincidence was observed; the only discrepancy occurred 
for the azimuth error in the front and back winds , where the simulated errors were 
smaller than the measured ones. This happened because of the assumption of the 
spatial coherence of the wind gust model, which generated a symmetric antenna 
load and, in consequence, a low azimuth servo error. This result indicates a need 
for upgrading the wind gust model to a spatially incoherent one that will reflect the 
real gusts in a more accurate manner. 

In order to design a controller with wind disturbance rejection properties, the 
wind disturbance should be known at the input to the antenna rate loop model. The 
second task, therefore, consists of developing a digital filter that simulates the wind 
gusts at the antenna rate input. This filter matches the spectrum of the measured 
servo errors. In this scenario, the wind gusts are generated by introducing white 
noise to the filter input. 


I. Introduction 

The steady-state wind pressure distribution on scaled antenna models was measured during wind 
tunnel experiments, 1,2,3 and their validity for actual field antennas was unknown. Wind field data 
collected recently at the DSS-13 antenna were used to evaluate the accuracy of the steady wind pressure 
measured in the wind tunnel [2]. A similar evaluation can be done for the time- varying part (gusts) of 
the wind. 

The wind gust analytical model, as developed in [1], is used to simulate the pointing errors of the DSN 
antennas. The model was developed using the wind tunnel data (as in Footnotes 1 through 3) and the 
Davenport spectra, but its accuracy was unverified. In this article, the wind measurements of servo errors 
obtained on January 24, 1994, at the DSS-13 antenna site, c.f. [2], were compared with the simulated 
servo errors. In most cases, the comparison shows satisfactory coincidence between the measured and the 
simulated data. 


1 N. L. Fox and B. Layman, Jr., “Preliminary Report on Paraboloidal Reflector Antenna Wind Tunnel Tests,” JPL Interoffice 
Memorandum CP-3 (internal document), Jet Propulsion Laboratory, Pasadena, California, 1962. 

2 N. L. Fox, “Load Distributions on the Surface of Paraboloidal Reflector Antennas,” JPL Interoffice Memorandum CP-4 
(internal document), Jet Propulsion Laboratory, Pasadena, California, 1962. 

* R. B. Blaylock, “Aerodynamic Coefficients for a Model of a Paraboloidal Reflector Directional Antenna Proposed for a 
JPL Advanced Antenna System,” JPL Interoffice Memorandum CP-6 (internal document), Jet Propulsion Laboratory, 
Pasadena, California, 1964. 
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Recently, the linear quadratic Gaussian (LQG) controller for the DSS-13 antenna was designed and 
tested (see [4]). This model-based controller used the identified DSS-13 antenna model based on field 
experiments [3]. This antenna model does not include the wind disturbances, which are necessary for the 
design of an improved LQG controller with wind disturbance rejection properties. For this purpose, the 
wind-measured data were used to create the wind disturbance input into the antenna rate-loop model and 
will serve as a base for the design of an improved controller with wind disturbance rejection properties. 


II. Evaluation of the Analytical Model 

In the field experiment, the servo errors due to wind gusts were measured for the elevation angles from 
11 to 89 deg and for the the yaw angles (antenna azimuth position with respect to the wind direction) 
from 0 to 360 deg. The servo errors from the analytical model are available for elevation angles of 60 and 
90 deg and for yaw angles of 0 (front wind), 90 (side wind), and 180 deg (back wind). The results are 
obtained in the form of the standard deviations of the measured servo error, typically of the length of 
8000 samples collected at a sampling time of 0.02 s. The results of field measurements and simulations 
are shown in Figs. 1 through 4, where “x” denotes field data and “o” denotes the analytical results. For 
the elevation servo error measurements, there were multiple collections of the field data for each elevation 
position. Thus, in this case, the maximal and minimal root-mean-square sums of the measured error are 
plotted with the gray area between them (see Figs. 1 through 3). For the azimuth errors, there was one 
collection of data, so the field errors do not include the gray area. 

The elevation servo error plots indicate that the analytical error lies within the gray area of the 
min-max measurements, while the results of the azimuth error show very close relationship between the 
measured and simulated standard deviations of the servo error for the side wind and a discrepancy for 
the front and back winds. In the latter case, the analysis underestimates the error because of the symmetry 



Fig. 1 . Standard deviation of the elevation encoder output due to 
40-km/h wind front gusts. 



Fig. 2. Standard deviation of the elevation encoder output due to 
40-km/h wind side gusts. 
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Fig. 3. Standard deviation of the elevation encoder output due to 
40-km/h wind back gusts. 



Fig. 4. Standard deviation of the elevation encoder output due to 
32-km/h wind gusts. 


and the coherence of the wind loading model. In reality, wind gusting is significantly unsymmetrical and 
is spatially uncorrelated. This discrepancy can be corrected by introducing an incoherent wind model 
using cross-spectra (see [5]). 


III. Wind Gust Model Derived From the Field Data 

LQG controllers developed for the DSN antennas (see [4]) are based on the antenna model obtained 
from the field testing rather than on its analytical model. In order to evaluate the controller wind distur- 
bance rejection properties as well as to improve these properties, one has to develop a wind disturbance 
model compatible with the antenna-identified model. 


The antenna rate-loop model was identified for the azimuth and elevation loops separately. The cross- 
coupling between the azimuth and elevation axes, and vice versa, was low and was, therefore, ignored. 
The input to the model is the rate command, and the output is the encoder reading. The rate command 
creates difficulties in implementation of the analytical wind gust model because the wind is modeled as 
pressure at the antenna structure and is not readily transformed into the rate command disturbance, but 
this can be done by using the measured servo errors due to the wind gust disturbance. 


Further, only the antenna model in azimuth is considered (the elevation model is developed similarly). 
The filter at the rate input will model the wind gusts (see Fig. 5). The white noise disturbance, w(t), 
of unit intensity at the input to the wind filter is assumed. The filter transfer function, G(s), is to be 
determined. The filter output, tc r (t), adds to the rate command, and it serves as the wind gust model. 
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RATE LOOP 
MODEL 


e(t) 


Fig. 5. Wind filter configuration. 


Let the servo error due to the wind gusts be e(t) and its spectrum be e(co). The servo error due to 
disturbance vi(t) is e s (t), and its spectrum is e s (u). The filter transfer function, G(uj), is determined such 
that the difference between the simulated and the measured power spectrum is minimized as follows: 

G(lj) such that || e(u>) - e s (u>) || is minimal (1) 

The filter that satisfies Condition (1) is called the wind filter. 

Let G r (u;) be the antenna rate-loop transfer function from the rate input to the encoder servo error. 
Then the simulated error due to wind gusts is obtained as 


e s (u;) = G t (u;)G(lj)w(ll>) 


( 2 ) 


In Eq. (2), the spectrum w(u) is constant (independent of frequency), and the transfer function G r (u>) 
is dominated by the antenna resonance frequencies. In this case, the magnitude of the filter transfer 
function can be assumed to be a smooth curve in the form of a shaped integrator, that is, 


k (Tis + l) 2 
“ s(T 2 s + l)(T 3 s+l) 2 


(3a) 


where the time constants are 


Ti 


1 

27r/i 


t 2 


1 

2 7T/ 2 


> 


T 3 


1 

2tt/ 3 


(3b) 


and /i = 2.2 Hz, / 2 = 7.0 Hz, and / 3 = 12.0 Hz ate the frequencies where the magnitude of the in- 
tegrator k/s is shaped. The frequencies /, and / 2 determine the bandwidth of resonance frequencies 
of the antenna, and the frequency f 3 is the cut-off frequency for the wind disturbances. In this trans- 
fer function, the only unknown parameter is the gain, k. The plot of G(u) is shown in Fig. 6 for k = 1. The 
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Fig. 6. Magnitude of the wind filter transfer function for k = 1 . 


choice of the transfer function shape as in Eq. (3) was done after the investigation of the more general 
case, where G(s) was a rational function of polynomials of order 5 or less. The performance errors for 
the polynomials were almost the same as for G(s) in Eq. (3). 


The wind model for the azimuth rate loop was determined for two antenna elevation positions: 60 and 
11 deg. For each elevation position, the wind from the front, side, and back was considered. The spectra 
of the azimuth encoder output, measured and simulated, are shown in Fig. 7 for an elevation angle of 
60 deg and a wind direction from the back of the antenna. The measured spectrum shows two resonances, 
at 1.7 and 4.2 Hz, and the spectrum from simulations has an additional resonance peak at 3.1 Hz. The 
spectra are coincidental at the first three frequencies. The time series of measured and simulated encoder 
outputs are shown in Figs. 8(a) and 8(b), respectively. They show the similarity; the difference between 
their standard deviations was less than 7 percent. The gain, fc, in this case was 0.0095. 



Similar results were obtained for other cases. Gain k for an 11-deg elevation angle was as follows: 


Front, wind Side wind Back wind 


0-0075 0.0079 0.0075 


Gain A; for a 60-deg elevation angle was as follows: 


34 





Front wind 

Side wind 

Back wind 

0.0095 

0.0096 

0.0092 


The tables show that for a given elevation angle the gains for front, back, and side winds are almost the 
same. Therefore, the wind filter is independent of wind direction; however, it depends on the antenna 
elevation angle. 


2i i i i r 
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Fig. 8. Azimuth encoder error due to 32-km/h wind gusts from the back at a 60-deg elevation 
angle: (a) measured and (b) simulated. 


IV. Conclusions 

The measured wind data at the DSS-13 site were used to verify the analytical wind model. The 
comparison showed that servo errors from the analytical model fall within the measured servo error 
boundaries. However, for the front and back winds, the simulated azimuth errors were smaller than the 
measured ones. This occurred because of the assumption of the spatial coherence of the wind gust model. 
The coherence caused a symmetric antenna load and, in consequence, a low azimuth servo error. This 
shortcoming indicates a need for upgrading the analytical wind model so that the spatially incoherent 
wind gust model reflects the real gusts in a more accurate manner. 

The measured wind data were also used to generate a new wind model more suitable for the design of 
a control system with wind disturbance rejection properties. The wind filter was obtained for the antenna 
azimuth model for different elevation angles and different wind directions such that the simulated servo 
error is close to the measured one. 
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Proper modeling of the Global Positioning System (GPS) satellite yaw attitude 
is important in high-precision applications. A new model for the GPS satellite yaw 
attitude is introduced that constitutes a significant improvement over the previously 
available model in terms of efficiency, flexibility, and portability. The model is 
described in detail , and implementation issues, including the proper estimation 
strategy, are addressed . The performance of the new model is analyzed, and an 
error budget is presented. This is the first self-contained description of the GPS 
yaw attitude model. 


I. Introduction 

On June 6, 1994, the U.S. Air Force implemented a yaw bias on most Global Positioning System 
(GPS) satellites. By January 1995, the implementation was extended to all the satellites except SVN 10. 
The yaw bias was introduced as a way to make modeling of the yaw attitude of the GPS satellites 
during shadow crossings possible [2]. The yaw attitude of a biased GPS satellite during eclipse seasons is 
markedly different from the yaw attitude of a noneclipsing satellite or from that of an unbiased satellite. 
The yaw attitude of the GPS satellite has a profound effect on precise applications. Mismodeling the 
satellite attitude can cause decimeter-level error in the positioning of ground stations with certain GPS- 
based techniques and skew media calibrations. This required the development of a special attitude model 
for biased GPS satellites. In addition to the yaw bias effects, that model also corrected other mismodeling 
that existed in the old model, namely, that of the “noon turn.” 

The first attitude model written for the biased constellation was made freely available to the GPS 
community in the form of a collection of FORTRAN routines [1 ] . For simplicity, this model is referred to 
in this article as GYM94 (for GPS Yaw Attitude Model — 94). GYM94 was implemented in JPL’s GIPSY 
software and, in various forms, in other high-precision geodetic packages. The model was successfully 
used within JPL’s routine processing of daily GPS orbits and ground station coordinates for the Inter- 
national Global Positioning System Service (IGS). The model had some drawbacks, though. Mainly, it 
was cumbersome to implement and very demanding of computer resources, namely, memory and central 
processing unit (CPU) time. 

In this article, we describe a new model for the GPS satellite attitude, referred to as GYM95. The 
model is analytic, in contrast to the numerical nature of GYM94, which required sequential processing 
in time. A time series of yaw rates estimated by the routine GPS processing at JPL will be analyzed to 
demonstrate the need to estimate the yaw rates. 
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II. Background 

The analysis that led to the implementation of the yaw bias on GPS satellites is described in Bar-Sever 
et al. [2]. A general description of the first yaw attitude model can also be found there. For completeness, 
we give here a brief summary. 

The nominal yaw attitude of a GPS satellite is determined by satisfying two constraints: first, that 
the navigation antennas point toward the geocenter and, second, that the normal to the solar array 
surface will be pointing at the Sun. To meet these two conditions, the satellite has to yaw constantly. 
The resulting yaw attitude algorithm is singular at two points — the intersections of the orbit with the 
Earth-Sun line. At these points, the yaw attitude is not single- valued, as any yaw angle allows optimal 
view of the Sun. In the vicinity of these singular points, the yaw rate of the spacecraft, required to 
keep track of the Sun, is unbounded. This singularity problem was largely ignored prior to the release of 
GYM94. While this mismodeling problem could be fixed easily through the realization of a finite limit on 
the spacecraft yaw rate, a bigger problem existed that could only be addressed by changing the attitude 
control subsystem (ACS) on board the spacecraft. The ACS determines the yaw attitude of the satellite 
by using a pair of solar sensors mounted on the solar panels. As long as the Sun is visible, the signal from 
the solar sensors is a true representation of the yaw error. During shadow, in the absence of sunlight, the 
output from the sensors is essentially zero and the ACS is driven in an open-loop mode by the noise in 
the system. It turns out that even a small amount of noise can be enough to trigger a yaw maneuver at 
maximum rate. To make it possible to model the yaw attitude of the GPS satellites, the ACS had to be 
biased by a small but fixed amount. Biasing the ACS means that the Sun sensor’s signal is superposed 
with another signal (the bias) equivalent to an observed yaw error of 0.5 deg (the smallest bias possible). 
As a result, during periods when the Sun is observed, the satellite yaw attitude will be about 0.5 deg in 
error with respect to the nominal orientation. During shadow, this bias dominates the open-loop noise 
and will yaw the satellite at full rate in the direction of the bias. Upon shadow exit, the yaw attitude of 
the satellite can be calculated, and the Sun recovery maneuver can also be modeled. 

GYM94 accounted for the yaw bias as well as the limit on the yaw rate. It computed the satellite yaw 
angle through numerical integration of a control law. Its output was a large file containing the yaw attitude 
history and, optionally, partial derivatives of the yaw attitude with respect to the yaw rate parameter. 
This file could later be interpolated to retrieve a yaw angle at the requested time. This process required 
relatively large amounts of computer memory and CPU time. In addition, the model’s complex control 
law — a simulation of the onboard attitude determination algorithm — did not allow much physical insight 
into the problem and was hard to tune. To overcome all these deficiencies, the GYM95 model was created. 
GYM94 was used in studies of GPS calibration for the DSN since September 1994, and the design of the 
new attitude model drew on the experience accumulated with GYM94. GYM95 is simple enough to be 
described by a small set of formulas, allowing easy implementation in different computing environments. 
Its analytic nature, as opposed to the numerical nature of GYM94, allows queries at arbitrary time points 
with great savings in computer resources. Finally, it allows more flexibility in tuning and adapting it to 
the changing conditions of the GPS constellation. 


III. The New Yaw Attitude Model (GYM95) 

A. Overview 

The yaw attitude of a GPS satellite can be divided into four regimes: nominal attitude, shadow 
crossing, postshadow maneuver, and noon turn. Most of the time (and for noneclipsing satellites all the 
time), the satellite is in the nominal attitude regime. The postshadow maneuver begins immediately after 
emerging from the Earth’s shadow and lasts until the satellite has regained its nominal attitude. This 
phase can last from 0 to 40 min. The noon-turn maneuver does not occur until the beta angle goes below 
about 5 deg and can last between 0 and 40 min. 
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We will start by defining a few important terms in Table 1 and the notation used, and then describe the 
yaw attitude during each of the four regimes, including the governing formulas. Finally, we will describe 
how to tie all the regimes together into one functional model and analyze any built-in errors. 


Table 1. Definition of terms. 


Term 

Definition 

Orbit midnight 

The point on the orbit furthest from the Sun. 

Orbit noon 

The point on the orbit closest to the Sun. 

Orbit normal 

The unit vector along the direction of the satellite’s angular momentum, treating 
the satellite as a point mass (equals position x velocity, where the order of 
the cross-product is important). 

Sun vector 

The direction from the spacecraft to the Sun. 

Beta angle 

The acute angle between the Sun vector and the orbit plane. It is defined as 
positive if the Sun vector forms an acute angle with the orbit normal and 
negative otherwise. 

Orbit angle 

The angle formed between the spacecraft position vector and orbit midnight, 
growing with the satellite’s motion. 

Yaw origin 

A unit vector that completes the spacecraft position vector to form an 
orthogonal basis for the orbit plane and is in the general direction of the 
spacecraft velocity vector. 

Spacecraft-fixed z-axis 

The direction of the GPS navigation antennas. 

Nominal spacecraft-fixed x-axis 

A unit vector orthogonal to the spacecraft- fixed z-axis 

such that it lies in the Earth-spacecraft-Sun plane and points in the general 
direction of the Sun (note that this definition is not single valued when the 
Earth, spacecraft, and Sun are collinear). 

Spacecraft-fixed x-axis 

A spacecraft-fixed vector, rotating with the spacecraft, such that far 
enough from orbit noon and orbit midnight, it coincides with the nominal 
spacecraft-fixed x-axis. Elsewhere, it is a rotation of the nominal 
spacecraft-fixed x-axis around the spacecraft-fixed z-axis. 

Nominal yaw angle 

The angle between the nominal spacecraft-fixed x-axis and the 

yaw-origin direction, restricted to be in [-180,180]. It is defined to have a sign 

opposite to that of the beta angle. 

Yaw angle 

The angle between the spacecraft-fixed x-axis and the yaw-origin direction, 
restricted to be in [—180,180], also termed “actual yaw angle.” 

Yaw error 

The difference between the yaw angle and the nominal yaw angle, restricted 
to be in [-180,180]. 

Midnight turn 

The yaw maneuver the spacecraft is conducting from shadow entry until it 
resumes nominal attitude sometime after shadow exit. 

Noon turn 

The yaw maneuver the spacecraft is conducting in the vicinity of orbit noon 
when the nominal yaw rate would be higher than the yaw rate the spacecraft is 
able to maintain. It ends when the spacecraft lesurnes nominal attitude. 

Spin-up/down time 

The time it takes for the spacecraft to spin up or down to its maximal 
yaw rate. The spacecraft is spinning down when it has to reverse its yaw rate. 


The notation used is as follows: 

/i = orbit angle 
(3 = beta angle 

E = Earth- spacecraft- Sun angle 
b = yaw bias inserted in the satellite ACS 
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B = actual yaw angle induced by b 
$ = actual yaw angle 
= nominal yaw angle 
t — current time, s 
ti = time of shadow entry 
t e = time of shadow exit 
t n — start time of the noon-turn maneuver 
t\ = spin-up /-down time 
\^ i = yaw angle upon shadow entry 
vj) e = yaw angle upon shadow exit 
R = maximal yaw rate of the satellite 
RR = maximal yaw-rate rate of the satellite 

Angle units, i.e., radians or degrees, will be implied by context. Radians will usually be used in 
formulas, and degrees will usually be used in the text. FORTRAN function names are used whenever 
possible with the implied FORTRAN functionality, e.g., ATAN2(a,b) is used to denote arc-tangent (a/b) 
with the usual FORTRAN sign convention. 

B. The Nominal Attitude Regime 

The realization of the two requirements for satellite orientation mentioned above yields the following 
formula for the nominal yaw angle: 

- ATAN2(-TAN(j3), SIN(/x)) + B(b, /3, p) (1) 


where (3 is the beta angle, /i is the orbit angle, measured from orbit midnight in the direction of motion, 
and B is the yaw bias (see below). It follows from this formula that the sign of the yaw angle is always 
opposite that of the beta angle. 

Ignoring the time variation of the slow-changing beta angle leads to the following formula for the 
yaw rate (there are simpler formulas, but they contain removable singularities that are undesirable for 
computer codes): 


= TANM « COSO.) x slNM ,; TA ri W + 6(‘. A ri < 2 > 

where fi varies little in time and can safely be replaced by 0.0083 deg/s. Notice that the sign of the 
nominal yaw rate is the same as the sign of the beta angle in the vicinity of orbit midnight (/x = 0). 

The singularity of these two formulas when (3 = 0 and fi — 0, 180 is genuine and cannot be removed. 

C. The Yaw Bias 

Like any medicine, the yaw bias has its side effects. Outside shadow, it introduces yaw errors that 
are actually larger that 0.5 deg. To fully understand this, we have to describe the ACS hardware, which 
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is beyond the scope of this article. The underlying reason is that the output of the solar sensor is 
proportional not to the yaw error but to its sine, and it is also proportional to the sine of the Earth- 
spacecraft-Sun angle, E. So, in order to offset a bias of b deg inserted in the ACS, the satellite has to 
actually yaw B deg, where B is given by: 


£(M,P 0 =B{b,E) = ASIN 


0.0175 x b 
SIN (E) 


( 3 ) 


The hardware-dependent proportionality factor is 0.0175, and the Earth-spacecraft-Sun angle, E , the 
beta angle, /3, and the orbit angle, /i, satisfy the following approximate relationship: 


COS(E) - COS(/?) x COS (/i) 


( 4 ) 


and E is restricted to [0,180]. Equation (3) becomes singular for E less then 0.5013 deg. This has no 
effect on the actual yaw because a small value of E implies that the spacecraft is in the middle of a 
midnight turn or a noon turn and is already yawing at full rate. The value of B does have a significant 
effect, though, on the timing of noon-turn entry and on the yaw angle shortly before that. For example, 
for E = 5 deg, which is the typical threshold value for noon-turn entry, the actual yaw bias is B ^ 6 deg. 

The bias rate, B , is given by 


B{b,(3,v) 


-0.0175 x b x COS(£) x COS {(3) x SIN(/x) x - QS(£?) - ^ SIN(E )3 


( 5 ) 


The ACS bias, 6, can be ±0.5 deg or 0 deg. With few exceptions, to be discussed below, the bias is always 
set to b = ~SIGN(0.5,/3) since this selection was found to expedite the Sun recovery time after shadow 
exit. 

D. The Shadow-Crossing Regime 

As soon as the Sun disappears from view, the yaw bias alone is steering the satellite. On most satellites, 
the yaw bias has a sign opposite to that of the beta angle. To correct for the bias-induced error, the 
satellite has to reverse its yaw rate upon shadow entry. For those satellites with bias of equal sign to 
that of the beta angle, there is no yaw reversal. The bias is large enough to cause the satellite to yaw 
at full rate until shadow exit, when the bias can be finally compensated. The yaw angle during shadow 
crossing depends, therefore, on three parameters: the yaw angle upon shadow entry, the yaw rate 
upon shadow entry, and the maximal yaw rate, R. Let t t be the time of shadow entry and let t be 
the current time, and define 


SIGN 

tl - SIGN {RR,b) 


to be the spin-up/-down time. Then the yaw angle during shadow crossing is given 


f \Iq ± 4 * i x (t — U) ± 0.5 x SIGN x (t — U 
\ ^ X tl + 0.5 x SIGN (RR, b) x t\ + SIGN (R, b) x {t - U - t x ) 


Using this formula, we avoid the singularity problem of the nominal attitude at midnight. 


by 


t < t{ ± t\ 
else 


( 7 ) 
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E. The Postshadow Maneuver 

This is the trickiest part of the yaw attitude model. The postshadow maneuver depends critically 
upon the yaw angle at shadow exit. The ACS is designed to reacquire the Sun in the fastest way possible. 
Upon shadow exit, the ACS has two options: One is to continue yawing at the same rate until the nominal 
attitude is resumed; the second is to reverse the yaw rate and yaw at full rate until the nominal attitude 
is resumed. In this model, we assume that the decision is based on the difference between the actual 
yaw angle and the nominal yaw angle upon shadow exit, and we denote this difference by D. If t e is the 
shadow-exit time, then 


D = ¥„(*«) - 9(t e ) - NINT ( *" ( * C 3 C o— x 360 


(8) 


and the yaw rate during the postshadow maneuver will be SIGN (R, D). 

Given the yaw angle upon shadow exit, the yaw rate upon shadow exit, SIGN (R, fo), and the yaw rate 
during the postshadow maneuver, we can compute the actual yaw angle during the postshadow maneuver 
by using Eq. (7) with the appropriate substitutions. This yields 

_ SIGN(i?, D) - SIGN(i?, b) 
tl ~ SIGN(fl#, D) 

. f + SIGN(i?, b) x (t — t e ) + 0.5 x SIGN (RR, D) x (t - t e ) 2 t<t e +ti 

\ + SIGN(J?, b) x + 0.5 X SIGN (RR, D)xt\ + SIGN(tf, D) x (t - t e - G) else 

( 10) 

The postshadow maneuver ends when the actual yaw attitude, derived from Eq. (10), becomes equal to 
the nominal yaw attitude. The time of this occurrence is computed in GYM95 by an iterative process 
that brackets the root of the equation $( t ) = ^(t), where the time dependence of ^(t) I s introduced 
by substituting fi = /x e + 0.0083 x (t — t c ) in Eq. (1). This equation can be solved as soon as the satellite 
emerges from shadow. Once the time of resuming nominal yaw is reached, the satellite switches back to 
that regime. 

F. The Noon-Turn Regime 

The noon-turn regime starts in the vicinity of orbit noon, when the nominal yaw rate reaches its 
maximal allowed value, and ends when the actual yaw attitude catches up with the nominal regime. 
First, we have to identify the starting point, and this can be done by finding the root, t n , of the equation 
'l'n(t) — -SIGN (i?., /3), where $ n (t) is the nominal yaw rate from Eq. (2). After the start of the noon 
turn, the yaw angle is governed by Eq. (7), again with the proper substitutions. This yields 


« = *n(tn)-SIGN(B,/3)x(t-t fl ) 


( 11 ) 


The end time is found by the same procedure that is used to find the end time of the postshadow 
maneuver. 

G. The Complete Model 

Satellite position and velocity, as well as the timing of shadow crossings, are required inputs to GYM95. 
The model is able to bootstrap, though, if these input values are unavailable far enough into the past. 
For example, if the satellite is potentially in the postshadow regime upon first query, there is a need to 
know the shadow-entry time so that all the inputs to Eqs. (9) and (10) be known. If this shadow-entry 
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time is missing from the input, the model can compute it approximately, as well as the shadow-exit time. 
Once all the timing information is available, yaw angle queries can be made at arbitrary time points. The 
model will decide the relevant yaw regime and compute the yaw angle using the correct formula. Given 
the above formulas, it is an easy matter to compute the partials of the yaw angle with respect to any 
parameter of the problem, the most important of which is the maximal yaw rate, R. 

H. Model Fidelity 

The fidelity of the model is a measure of how accurately it describes the true behavior of the satellite. 
This is hard to measure because there is no high-quality telemetry from the satellite and because the 
estimated value of the main model parameter, namely, the yaw rate, depends on many other factors beside 
the attitude model itself: data, estimation strategy, and other models for the orbit and the radiometric 
measurements. Nevertheless, based on the experience accumulated thus far with this model and its 
predecessor, GYM94, it is possible to come up with an educated guess of the inaccuracy of GYM95. 

The nominal attitude regime is believed to be very accurate. The only source of error is mispointing 
of the satellite, which is poorly understood and relatively small (of the order of 1 deg around the pitch, 
yaw, and roll axes). Compensations for the dynamic effect of this error source are discussed in [3] and 
[4], where they are treated, properly, within the context of the solar pressure model. 

Modeling the midnight turn accurately is difficult. Inherent uncertainties like the exact shadow-entry 
and -exit times are a constant error source. Inaccuracies in shadow-entry time are more important than 
inaccuracies in shadow-exit time because errors in the former are propagated by the model throughout the 
midnight-turn maneuver. In contrast, error in the shadow-exit time will affect the postshadow maneuver 
only. Either way, the inaccuracy will be manifested through a constant error in the yaw angle, something 
that can be partially compensated through the estimation of the yaw rate. The length of the penumbra 
region is usually about 60 s. Sometime during this period, the yaw bias kicks in. GYM95 puts that time 
midway into penumbra. The maximum timing error is, therefore, less than 30 s. A worst-case scenario, 
ignoring the short spin-up/-down period and using a yaw rate of 0.13 deg/s, will give rise to a constant 
yaw error of 30 x 0.13 « 4 deg throughout the midnight turn. A more realistic estimate is 3 deg, even 
before applying yaw rate compensation, after which the rms error will remain the same, but the mean is 
expected to vanish. Another error source is the uncertainty in the value of the maximal yaw-rate rate, 
RR. This parameter is weakly observable and, therefore, hard to estimate. The nominal value used in 
GYM95 is 0.00165 deg/s 2 for Block II A satellites and 0.0018 deg/s 2 for Block II satellites, and it should 
be at least 70-percent accurate. The long-term effects of a yaw-rate rate error can be computed from the 
second part of Eq. (7) as 


*(RR) = 


x SIGNER, b) - 0.5 x^- 0.5 x SIGN( fl, b) 2 
sign (RR, b) 


A worst-case scenario assuming = — STGN(i?., b ) =0.13 and 30- percent error in the yaw-rate rate would 
give rise to a yaw error of about 5 deg. These assumptions also imply a very short shadow duration so 
the error will not be long lasting. For long shadow events, w 0, and the resulting yaw error is about 
1 deg. Again, this error can be partially offset by estimating the yaw rate. 

The main error source for the noon turn is the timing uncertainty of the onset of the maneuver. This 
uncertainty is not expected to be larger than 2 min. A 2-min error will cause a constant yaw error of about 
15 deg, assuming a yaw rate of 0.13 deg/s. The relatively short duration of the noon turn diminishes 
somewhat the effects of such a large error. Estimating the yaw rate will decrease the error further. 

The value of the yaw rate is not considered here as an error source. Any nominal value stands to be 
at least 10 percent in error (see below). Since errors due to yaw rate grow in time, this parameter must 
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be estimated or, alternatively, a previously estimated value should be used. For example, an error of 
0.01 deg/s in the yaw rate will give rise to a 30-deg error in yaw at the end of a 50-min shadow event. 

Although unlikely, errors from different sources can add up. In that case, the maximal error for each 
regime is as follows: 2 deg for the nominal yaw regime, 9 deg for the midnight-turn regime, and 15 deg 
for the noon-turn regime. Typical errors are expected to be less than half of these values. 


IV. The Estimated Yaw Rates 

As part of the implementation of the GYM models at JPL, the yaw rates of all eclipsing satellites 
are estimated for every midnight turn and every noon turn. In JPL’s GIPSY software, this is done 
by treating the yaw rate as a piece-wise constant parameter for each satellite. The parameter value is 
allowed to change twice per revolution, midway between noon and midnight. Since a small error in the 
yaw rate can cause a large yaw error over time and, since our a priori knowledge of the yaw rate is not 
sufficiently accurate, we found it necessary to iterate on the yaw rate value. JPL routinely publishes the 
final estimates for the yaw rates as daily text files. Unfortunately, due to a software bug, the archived yaw 
rates for dates prior to February 16, 1995, were in error. This leaves a period of about 2 months when the 
estimated yaw rates are available. Figure 1 depicts the estimated yaw rates for each eclipsing satellite, 
for each midnight turn, and for each noon turn, from February 16 to April 26, 1995. The accuracy 
of the estimates depends on the amount of data available during each maneuver and this, in turn, is 
proportional to the duration of the maneuver. The longer the maneuver, the better the estimate. The 
effect of a reduced estimation accuracy during short maneuvers is mitigated by the fact that the resulting 
yaw error is also proportional to the duration of the maneuver. For long maneuvers, e.g., a midnight turn 
at the middle of the eclipse season, the estimates are good to 0.002 deg/s, which leads to a maximal yaw 
error of about 6 deg. A similar error level is expected for short maneuvers. Noon turns occur only during 
the middle part of the eclipse season. In Fig. 1, they can be distinguished from midnight-turn rates by the 
larger formal error associated with them, since they are typically short events of 15- to 30- min duration. 
As a result, the scatter of the noon-turn rates is larger than that of the midnight-turn rates. Toward 
the edges of the eclipse season, the quality of the yaw rate estimates drops, again because of the short 
duration of the shadow events. The most striking feature in Fig. 1 is the discontinuity of the estimated 
yaw rates in the middle of the eclipse season, corresponding to the beta angle crossing zero. No plausible 
explanation is currently available for this jump. SVN 29 is the only satellite that does not have a jump 
discontinuity; this is also the only satellite that does not undergo a bias switch in the middle of the eclipse 
season. SVN 31 is the only satellite with a jump from high yaw rates to low yaw rates as the beta angle 
transitions from positive to negative. There is nothing otherwise special about SVN 31. The ratio of the 
high yaw-rate values to the low yaw-rate values is about 1.3 for all satellites. 

Within each half of the eclipse season, the midnight yaw rates are fairly constant, varying by 10 percent 
or less. The noon-turn yaw rates seem to be more variable. This is not only a consequence of the weak 
observability but also of the fact that the spacecraft is subject to a varying level of external torque during 
the noon turn as the eclipse season progresses. 

The modeling of the postshadow maneuver is a problem for which a satisfactory solution has not 
yet been found. The source of the problem is that the presence of the postshadow regime makes the 
estimation of the yaw rate into a nonlinear problem. There is always a critical value of the yaw rate 
such that, for higher values, the spacecraft will reverse its yaw upon shadow exit and, tor lower values, 
the spacecraft will retain its yaw rate until the end of the midnight turn. If this critical value falls in 
the range of feasible yaw rates— -which it often does— it becomes very hard to figure out what kind of 
maneuver the satellite is doing upon shadow exit. To avoid this postshadow ambiguity, we have been 
rejecting measurement data from shadow exit until about 30 min thereafter. 
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Fig. 1. Estimated yaw rates with their formal errors versus GPS week for coplanar (C-plane) (a) SVN 28, 
(b) SVN 31, (c) SVN 36, and (d) SVN 37 and for F-plane (e) SVN 18, (f) SVN 26, (g) SVN 29, and (h) SVN 32. 
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We describe a novel oscillator that converts continuous light energy into sta- 
ble and spectrally pure microwave signals. This light-induced microwave oscillator 
(LIMO) consists of a pump laser and a feedback circuit, including an intensity 
modulator, an optical fiber delay line, a photodetector, an amplifier, and a filter . 
We develop a quasilinear theory and obtain expressions for the threshold condition, 
the amplitude, the frequency ; the line width, and the spectral power density of 
the oscillation. We also present experimental data to compare with the theoret- 
ical results. Our findings indicate that the LIMO can generate ultrastable, spec- 
trally pure microwave reference signals up to 75 GHz with a phase noise lower than 
— 140 dBc/Hz at 10 kHz. 


I. Introduction 

Oscillators are devices that convert energy from a continuous source to a periodically varying signal. 
They represent the physical realization of a fundamental basis of all physics, the harmonic oscillator, and 
they are perhaps the most widely used devices in modern day society Today a variety of oscillators — 
mechanical [1] (such as the pendulum), electromagnetic (such as LC circuit [2,3] and cavity-based [4]), and 
atomic (such as masers [5] and lasers [6]) — provides a diverse range in the approximation to the realization 
of the ideal harmonic oscillator. The degree of the spectral purity and stability of the output signal of the 
oscillator is the measure of the accuracy of this approximation and is fundamentally dependent on the 
energy storage ability of the oscillator, determined by the resistive loss (generally frequency dependent) 
of the various elements in the oscillator. 


An important type of oscillator widely used today is the electronic oscillator. The first such oscillator 
was invented by L. De Forest [2] in 1912, shortly after the development of the vacuum tube. In this 
triode-based device known as the van der Pol oscillator [3], the flux of electrons emitted by the cathode 
and flowing to the anode is modulated by the potential on the intervening grid. This potential is derived 
from the feedback of the current in the anode circuit containing an energy storage element (i.e., the 
frequency-selecting LC filter) to the grid, as shown in Fig. 1(a). Today the solid-state counterparts 
of these valve oscillators based on transistors are pervasive in virtually every application of electronic 
devices, instruments, and systems. Despite their widespread use, electronic oscillators, whether of the 
vacuum-tube or the solid-state variety, are relatively noisy and lack adequate stability for applications 
where very high stability and spectral purity are required. The limitation to the performance of electronic 
oscillators is due to ohmic and dispersive losses in various elements in the oscillator, including the LC 
resonant circuit. 
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(a) 


/ PLATE CURRENT 




Fig. 1. Comparison of two types of oscillators: (a) van der Pol oscillator and (b) light-induced microwave oscillator. 


For approximately the past 50 years, the practice of reducing the noise in the electronic oscillator by 
combining it with a high-quality factor (Q) resonator has been followed to achieve improved stability and 
spectral purity. The Q is a figure of merit for the resonator given by Q — 2n fr^ where is the energy 
decay time that measures the energy storage ability of the resonator and / is the resonant frequency. 
High-Q resonators used for stabilization of the electronic oscillator include mechanical resonators, such 
as quartz crystals [7,8]; electromagnetic resonators, such as dielectric cavities [9]; and acoustic [10] and 
electrical delay lines, where the delay time is equivalent to the energy decay time, r^, and determines 
the achievable Q. This combination with a resonator results in hybrid-type oscillators referred to as 
electromechanical, electromagnetic, or electro-acoustic, depending on the particular resonator used with 
the oscillator circuit. The choice of the particular resonator is generally determined by a variety of factors, 
but for the highest achievable Q’s at room temperatures, the crystal quartz is the resonator of choice for 
the stabilization of the electronic oscillator. However, because quartz resonators have only a few high-Q 
resonant modes at low frequencies [7,8], they have a limited range of frequency tunability and cannot be 
used to directly generate high-frequency signals. 


In this article, we introduce a novel photonic oscillator [11,12] characterized by spectral purity and 
frequency stability rivaling the best crystal oscillators. This oscillator, shown schematically in Fig. 1(b), 
is based on converting the continuous light energy from a pump laser to radio frequency (RF) and 
microwave signals, and thus we refer to it with the acronym LIMO, for the light-induced microwave 
oscillator. The LIMO is fundamentally similar to the van der Pol oscillator, with photons replacing 
the function of electrons, an electro-optic (E/O) modulator replacing the function of the grid, and a 
photodetector replacing the function of the anode. The energy storage function of the LC circuit in the 
van der Pol oscillator is replaced with a long fiber-optic delay line in the LIMO. 

Despite this close similarity, the LIMO is characterized by significantly lower noise and very high 
stability, as well as by other functional characteristics that are not achieved with the electronic oscillator. 
The superior performance of the LIMO results from the use of electro-optic and photonic components that 
are generally characterized by high efficiency, high speed, and low dispersion in the microwave frequency 
regime. Specifically, currently there are photodetectors available with quantum efficiency as high as 
DO percent that can respond to signals with frequencies as high as 110 GHz [13]. Similarly, E/O modulators 
with a 75-GHz frequency response are also available [14]. Finally, the commercially available optical fiber, 
which has a small loss of 0.2 dB/km for 1550 nm light, allows long storage time of the optical energy 
with negligible dispersive loss (loss dependent on frequency) for the intensity modulations at microwave 
frequencies. 
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The LIMO may also be considered as a hybrid oscillator in so far as its operation involves both light 
energy and microwave signals. Nevertheless, as a hybrid oscillator, the LIMO is unique in that its output 
may be obtained both directly as a microwave signal or as intensity modulation of a optical carrier. This 
property of the LIMO is quite important for applications involving optical elements, devices, or systems 

in]. 

The ring configuration, consisting of an electro-optic modulator that is fed back with a signal from the 
detected light at its output, has been previously studied by a number of investigators interested in the 
nonlinear dynamics of bistable optical devices [15-19]. The use of this configuration as a possible oscillator 
was first suggested by A. Neyer and E. Voges [20]. Their investigations, however, were primarily focused on 
the nonlinear regime and the chaotic dynamics of the oscillator. This same interest persisted in the work of 
T. Aida and P. Davis [21], who used a fiber wave guide as a delay line in the loop. Our studies, by contrast, 
are specifically focused on the stable oscillation dynamics and the noise properties of the oscillator. The 
sustainable quasilinear dynamics, both in our theoretical and experimental demonstrations, are arrived 
at by the inclusion of a filter in the feedback loop to eliminate harmonics generated by the nonlinear 
response of the E/O modulator. This approach yields stable, low-noise oscillations that closely support 
the analytical formulation presented here. 

In this article, we first describe the oscillator and identify the physical basis for its operation. We then 
develop a quasilinear theory for the oscillator dynamics and the oscillator noise. Results of the theory 
are then compared with experimental results. 


il. Description of the Oscillator 

The LIMO utilizes the transmission characteristics of a modulator together with a fiber-optic delay 
line to convert light energy into stable, spectrally pure RF /microwave reference signals. A detailed view 
of the construction of the oscillator is shown schematically in Fig. 2. In this depiction, light from a laser 
is introduced into an E/O modulator, the output of which is passed through a long optical fiber, and 
detected with a photodetector. The output of the photodetector is amplified and filtered and fed back to 
the electric port of the modulator. This configuration supports self-sustained oscillations at a frequency 
determined by the fiber delay length, bias setting of the modulator, and the bandpass characteristics of 
the filter. It also provides for both electrical and optical outputs, a feature which would be of considerable 
advantage to photonics applications. 

We use a regenerative feedback model to analyze the spectral properties of the LIMO. Similar methods 
have been successfully used to analyze lasers [5] and surface acoustic wave oscillators [22]. The conditions 
for self sustained oscillations include coherent addition of partial waves each way around the loop and a 
loop gain exceeding losses for the circulating waves in the loop. The first condition implies that all signals 
that differ in phase by some multiple of 2n from the fundamental signal may be sustained. Thus, the 
oscillation frequency is limited only by the characteristic frequency response of the modulator and the 
setting of the filter, which eliminates all other sustainable oscillations. The second condition implies that 
with adequate light input power, self sustained oscillations may be obtained, without the need for the 
RF/microwave amplifier in the loop. These characteristics, which are expected based on the qualitative 
analysis of the oscillator dynamics, are mathematically derived in the following sections. 


III. Quasilinear Theory of the LIMO 

In the following sections, we introduce a quasilinear theory to study the dynamics and noise of a 
LIMO. In the discussion, we assume that the E/O modulator in the oscillator is of the Mach-Zehnder 
type. However, the analysis of oscillators made with different E/O modulators may follow the same 
procedure. The flow of the theory is as follows: First, the open-loop characteristics of a photonic link 
consisting of a laser, a modulator, a fiber delay, and a photodetector are determined. We then close the 
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Fig. 2. Detailed construction of a LIMO. Optical injection and RF injection ports are 
supplied for synchronizing the oscillator with an external reference by either optical 
injection locking or electrical injection locking [12]. The bias port and the fiber stretcher 
can be used to fine tune the oscillation frequency [12]. Noise in the oscillator can be 
viewed as being injected from the input of the amplifier. 


loop back into the modulator and invoke a quasilinear analysis by including a filter in the loop. This 
approach leads us to a formulation for the amplitude and the frequency of the oscillation. In the next 
step, we consider the influence of the noise in the oscillator, again assisted by the presence of the filter, 
which limits the number of circulating Fourier components. We finally arrive at an expression for the 
spectral density of the LIMO that would be suitable for experimental investigations. 

A. Oscillation Threshold 

The optical power from the E/O modulator's output port that forms the loop is related to an applied 
voltage V in (t) by 


P(t) = 



Tf sill 7 T 


V<n(0 

K 



( 1 ) 


where a is the fractional insertion loss of the modulator, V n is its half-wave voltage, Vb is its bias voltage, 
P 0 is the input optical power, and 77 determines the extinction ratio of the modulator by (1 4- 77) /(I - 77). 

If the optical signal P(t) is converted to an electric signal by a photodetector, the output electric signal 
after an RF amplifier is 


Vout(t) = pP(t)RG A = V p h{ 1 - 77 sin 7r 


Vin(t) Vb 

Vk V ; 


( 2 ) 


where p is the responsivity of the detector, R is the load impedance of the photodetector, G A is the 
amplifier’s voltage gain, and V ph is the photovoltage, defined as 


v ph = RG a = I vh RG A 


( 3 ) 


with I ph = aP 0 p /2 as the photocurrent. The LIMO is formed by feeding the signal of Eq. (2) back to the 
RF input port of the E/O modulator. Therefore, the small signal open- loop gain, Gs, of the LIMO is 
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V; 7l =0 


( 4 ) 


G s = 


dVput 

~dV~ 



The highest small signal gain is obtained when the modulator is biased at quadrature, that is, when 
Vb = 0 or Vn- From Eq. (4), one may see that Gs can be either positive or negative, depending on 
the bias voltage. The modulator is said to be positively biased if Gs > 0; otherwise it is negatively 
biased. Therefore, when V B = 0, the modulator is biased at negative quadrature, while when V B = V n , 
the modulator is biased at positive quadrature. Note that in most externally modulated photonic links 
the E/O modulators can be biased at either positive or negative quadrature without affecting their 
performance. However, as will be seen next, the biasing polarity will have an important effect on the 
operation of the LIMO. 

In order for the LIMO to oscillate, the magnitude of the small signal open-loop gain must be larger 
than unity. From Eq. (4), we immediately obtain the oscillation threshold of the LIMO to be 

Kr 

Vph 7TT7|cos(?rV B /%)| ^ 

For the ideal case in which 77 = 1 and V B — 0 or V n , Eq. (5) becomes 

Vph = — (6) 

7 r 

It is important to notice from Eqs. (3) and (6) that the amplifier in the loop is not a necessary condition 
for oscillation. So long as I p hR > / 7r is satisfied, no amplifier is needed ( Ga = 1)- It is the optical 

power from the pump laser that actually supplies the necessary energy for the photonic oscillator. This 
property is of practical significance because it enables the LIMO to be powered remotely using an optical 
fiber. Perhaps more significantly, however, the elimination of the amplifier in the loop also eliminates 
the amplifier noise, resulting in a more stable oscillator. For a modulator with a V n of 3.14 V and 
an impedance, R , of 50 a photocurrent of 20 mA is required for sustaining the photonic oscillation 
without an amplifier. This corresponds to an optical power of 25 mW, assuming the responsivity, p, of 
the photodetector to be 0.8 A/W. 

B. Linearization Of the E/O Modulator’s Response Function 

In general, Eq. (2) is nonlinear. If the electrical input signal Vi n {t) to the modulator is a sinusoidal 
wave with an angular frequency of u, an amplitude of V 0 , and an initial phase of /?, 


V in (t) = V 0 sm(ust-h/3) (7) 

then the output at the photodetector, V out (t ) , can be obtained by substituting Eq. (7) in Eq. (2) and 
expanding the left-hand side of Eq. (2) with Bessel functions: 


Fout(^) — Vph 



cos(2 mut 4- 2m/?) 


- 2r] cos 




sin [(2m + 1 )ujt + (2m + 1)/?] 


( 8 ) 
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It is clear from Eq. (8) that the output contains many harmonic components of lj. 

The output can be linearized if it passes through an RF filter with a bandwidth sufficiently narrow to 
block all harmonic components. The linearized output can be obtained easily from Eq. (8): 


Vout(t) = G(V 0 )V in (t) 


where the voltage gain coefficient, G(Vq), is defined as 


2V 

G(V 0 ) = G s —^J l 

7T Vq 



(9) 


(10a) 


It can be seen that the voltage gain, G(V 0 ), is a nonlinear function of the input amplitude, V Q , and its 
magnitude decreases monotonically with V Q . However, for a small enough input signal (V Q << V n and 
J^TiVo/V-n) - 7vV 0 /2V n ), we can recover from Eq. (10a) the small signal gain: G(V 0 ) — Gs- 

If we expand the left-hand side of Eq. (2) with Taylor series, the gain coefficient can be obtained as 


G(V 0 ) = G s 



(ttVo 

\2Vn 




(10b) 


It should be kept in mind that, in general, G ( V () ) is also a function of the frequency, u, of the input 
signal, because V p h is linearly proportional to the gain of the RF amplifier and the responsivity of the 
photodetector, which are all frequency dependent. In addition, the V n of the modulator is also a function 
of the input RF frequency. Furthermore, the frequency response of the RF filter in the loop can also be 
lumped into G{V 0 ). In the discussions below, we will introduce a unitless complex filter function, F(cj), 
to explicitly account for the combined effect of all frequency-dependent components in the loop while 
treating G(V 0 ) to be frequency independent: 

F{u) = F(u>)e t<t>{u] (11) 

where 0(u;) is the frequency- dependent phase caused by the dispersive component in the loop and F(&) 
is the real normalized transmission function. Now Eq. (9) can be rewritten in complex form as 


V out (t) ^ F(u;)G(V 0 )V in (^t) 


(12) 


where V in and V out are complex input and output voltages. Note that although Eq. (12) is linear, the 
nonlinear effect of the modulator is not lost — it is contained in the nonlinear gain coefficient, G(V a ). 


C. Oscillation Frequency and Amplitude 

In this section, we derive the expressions for the amplitude and frequency of the LIMO. Like other 
oscillators, the oscillation of a LIMO starts from noise transient, which is then built up and sustained 
with feedback at the level of the oscillator output signal. We derive the amplitude of the oscillating 
signal by considering this process mathematically. The noise transient can be viewed as a collection of 
sine waves with random phases and amplitudes. To simplify our derivation, we use this noise input with 
the linearized expression of Eq. (12) for the loop response. Because Eq. (12) is linear, the superposition 
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principle holds, and we can analyze the response of the LIMO by first inspecting the influence of a 
single-frequency component of the noise spectrum: 


V in (w,t) = V in (ui)e iuJt 


(13) 


where V in (u)) is a complex amplitude of the frequency component. 

Once the noise component of Eq. (13) is in the oscillator, it would circulate in the loop, and the 
recurrence relation of the fields from Eq. (12) is 


Vnfat) = F(u)G{V 0 )V n ^t^ r') 


(14) 


where r' is the time delay resulting from the physical length of the feedback and n is the number of times 
the field has circulated around the loop, with V n= o(cj,t) = V in (u>,t). In Eq. (14), the argument V Q in 
G(V 0 ) is the amplitude of the total field (the sum of all circulating fields) in the loop. 

The total field at any instant of time is the summation of all circulating fields. Therefore, with the 
input of Eq. (13) injected in the oscillator, the signal measured at the RF input to the modulator for the 
case when the open-loop gain is less than unity can be expressed as 


V(w,t) = G a V in (u>) Y, [F(u)G(V 0 )] n e iu{t - nT,) 

71 = 0 


GgVin^Y 

1 - F(u)G(V 0 )e-^ T ' 


(15) 


For a loop gain below threshold and with a small V a , G(V 0 ) is essentially the small signal gain, Gs , given 
by Eq. (4). 

The corresponding RF power of the circulating noise at frequency to is, therefore, 

|VW)| a <%|VinMI 2 /(2fl) , 1(n 

(u} 2 R 1 + \F(w)G(V 0 )\ 2 - 2F’(w)|G(C,)| cos[u;r' + <t>(u>) + 4>o}' K 1 


where 0 O = 0 if G(V 0 ) > 0 and (j) 0 = r if G{V 0 ) < 0. 

For a constant V in ( cj), the frequency response of a LIMO has equally spaced peaks similar to those of 
a Fabry-Perot resonator, as shown in Fig. 3. These peaks are located at the frequencies determined by 


u^t' + (p{^k) + 4>o ~ 2/c7t A; = 0, 1,2,--- (17) 

where k is the mode number. In Fig. 3, each peak corresponds to a frequency component resulting from 
the coherent summation of all circulating fields in the loop at that frequency. As the open-loop gain 
increases, the magnitude of each peak becomes larger and its shape becomes sharper. These peaks are 
the possible oscillation modes of the LIMO. When the open- loop gain is larger than unity, each time a 
noise component at a peak frequency travels around the loop, it is amplified and its amplitude increases 
geometrically — an oscillation is started from noise. 

Because an RF filter is placed in the loop, the gain of only one mode is allowed to be larger than unity, 
thus selecting the mode that is allowed to oscillate. Because of the nonlinearity of the E/O modulator 
or the RF amplifier, the amplitude of the oscillation mode cannot increase indefinitely. As the amplitude 
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Fig. 3. Illustration of the oscillator's output spectra below and above the threshold. 


increases, higher harmonics of the oscillation will be generated by the nonlinear effect of the modulator or 
the amplifier, at the expense of the oscillation power, and these higher harmonics will be filtered out by 
the RF filter. Effectively, the gain of the oscillation mode is decreased according to Eq. (10) until the gain 
is, for all practical measures, equal to unity, and the oscillation is stable. As will be shown later, because 
of the continuous presence of noise, the closed-loop gain of an oscillating mode is actually less than unity 
by a tiny amount on the order of 10 — 10 , which ensures that the summation in Eq. (15) converges. 

In the discussion that follows, only one mode k is allowed to oscillate, and so we will denote the 
oscillation frequency of this mode as f osc or uj 0 sc{u 0 sc — 27 r/ 05C ), its oscillation amplitude as V osc1 and its 
oscillation power as P 0 sc(Posc = Vosc! 2ft)- In this case, the amplitude, Vo, of the total field in Eq. (16) 
is just the oscillation amplitude, V osc , of the oscillating mode. If we choose the transmission peak of the 
filter to be at the oscillation frequency, cj 05C , and so F(lj osc ) — 1, the oscillation amplitude can be solved 
by setting the gain coefficient, |G(\^, SC )|, in Eq. (16) at unity. From Eq. (10a), this leads to 


J i 



1 7 tV osc 

2\G S \ V n 


(18a) 


In deriving Eq. (18a), we have assumed that the RF amplifier in the loop is linear enough that the 
oscillation power is limited by the nonlinearity of the E/O modulator. The amplitude of the oscillation 
can be obtained by solving Eq. (18a) graphically, and the result is shown in Fig. 4(a). Note that this 
result is the same as that obtained by Neyer and Voges [20] using a more complicated approach. 

If we use Eq, (10b), we can obtain the approximated solutions of the oscillation amplitude: 


Vn, 




i 

W\ 


thirdorder expansion (18b) 


v oat 


2V3K ^ 


1 

7s 



1/2 


fifthorder expansion 


(18c) 


The threshold condition of |Gs| > 1 is clearly indicated in Eqs. (18b) and (18c). Figure 4(a) shows 
the normalized oscillation amplitude as a function of |Gs|, obtained from Eqs. (18a), (18b) and (18c), 
respectively. Comparing the three theoretical curves, one can see that, for |Gs| < 1.5, the third-order 
expansion result is a good approximation. For |G$| < 3, the fifth-order expansion result is a good 
approximation. 
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Fig. 4. The normalized oscillation amplitude of a LIMO as a function 
of small signal gainJG^: (a) theoretical calculation using Eqs. (18a), 
(18b), and (18c) and (b) experimental data and curve fitting to Eqs. 
(18b) and (18c). 


The corresponding oscillation frequency, f osc = = Uk/ 2tt, can be obtained from Eq. (17) as 


_ , (* + 1/2) 
- Jk ~ 

T 

for G(V 0SC ) < 0 

(19a) 

= fk=- 

T 

for G(V r osc ) > 0 

(19b) 


where r is the total group delay of the loop, including the physical length delay, r', of the loop and the 
group delay resulting from dispersive components (such as an amplifier) in the loop, and it is given by 


T = T + 


d(f)(uj) 


duj 


( 20 ) 


For all practical purposes, J\{^V osc IV^) > 0 or V osc /Vn < 1.21, and the sign of G(V osc ) is determined by 
the small signal gain, G s . It is interesting to notice from Eqs. (19a) and (19b) that the oscillation frequency 
depends on the biasing polarity of the modulator. For negative biasing (G s < 0), the fundamental 
frequency is l/(2r), while for positive biasing (G s > 0), the fundamental frequency is doubled to 1/r. 
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D. The Spectrum 

The fundamental noise in a LIMO consists of the thermal noise, the shot noise, and the laser’s intensity 
noise, which for the purpose of analysis can be viewed as all originating from the photodetector. Since 
the photodetector is directly connected to the amplifier, the noise can be viewed as entering the oscillator 
at the input of the amplifier, as shown in Fig. 2. 

We compute the spectrum of the oscillator signal by determining the power spectral density of noise 
in the oscillator. Let Pn{u)) be the power density of the input noise at frequency <j; we have 


p N (uj)Af = 


lUnMI 2 

27? 


( 21 ) 


where Af is the frequency bandwidth. Substituting Eq. (21) in Eq. (16) and letting F(lj osc ) = 1, we 
obtain the power spectral density of the oscillating mode k: 


SrfU') 


P(f ) 

A fPosc 


PnG'a/ P osc 

1 + \F(f')G(V osc )\i - 2F(f')\G(V osc )\cos(27Tf'T) 


( 22 ) 


where = — uj osc ) /2tx is the frequency offset from the oscillation peak, / 05C « In deriving Eq. (22), 

both Eqs. (17) and (20) are used. 

By using the normalization condition 


/ oo /*l/2r 

SrfWW « / S RF (f')df = 1 (23) 

-oo «/— l/2r 

we obtain 

1 - |G(K osc )| 2 « 2 [1 - |G(V osc )|] = (24) 

Note that in Eq. (23) we have assumed that the spectral width of the oscillating mode is much smaller 
than the mode spacing, 1/r, of the oscillator, so that the integration over 1/r is sufficiently accurate. In 
addition, in the derivation, we have assumed that |F(/')| ~ 1 in the frequency band of integration. 

Typically, pjy ~ 10“ 17 mW/Hz, P osc ~ 10 mW, G\ ~ 100, and r ~ 10 -G s. From Eq. (24), one can 
see that the closed-loop gain, IG^V^c)!, of the oscillating mode is less than unity by an amount of 10“ 10 . 
Therefore, the equation 1(2(1^^) | — 1 is sufficiently accurate for calculating the oscillation amplitude, 
Vosc , as in Eqs. (18a), (18b), and (18c). 

Finally, substituting Ecp (24) in Eq. (22), we obtain the RF spectral density of the LIMO: 


SiifU') 


6 

(2 - 6/r) - 2yJ\ - 6/t cos(2nf'T) 


(25) 


where 6 is defined as 
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( 26 ) 


3 _ PnG\ 

Pose 

As mentioned before, p^ is the equivalent input noise density injected into the oscillator from the input 
port of the amplifier, and P osc jG 2 A is the total oscillating power measured before the amplifier. Therefore, 
6 is the input noise- to-signal ratio to the oscillator. 

For the case where 2 7t/'t << 1, we can simplify Eq. (25) by expanding the cosine function in Taylor 
series: 


SrfU') 


6 

(< 5 / 2 t ) 2 + ( 2 t ) 2 ( t /') 2 


(27) 


Equation (27) is a good approximation even for 27 r/'r = 0.7, at which value the error resulting from 
neglecting the higher-order terms in Taylor expansion is less than 1 percent. It can be seen from Eq. (27) 
that the spectral density of the oscillating mode is a Lorentzian function of frequency. Its full width at 
half maximum (FWHM), A /fwhm, is 


A Ifwhm = 


27 r r 2 


_L G aPn 

27 r t 2 P osc 


(28a) 


It is evident from Eq. (28a) that A Jfwhm is inversely proportional to the square of the loop delay time 
and linearly proportional to the input noise- to- signal ratio, 6. For a typical <5 of 10 -16 /Hz and a loop 
delay of 100 ns (20 m), the resulting spectral width is submilli-Hertz. The fractional power contained in 
A Ifwhm is A /fwhmS'/jf(0) = $4 percent. 

From Eq. (28a), one can see that for fixed px and Ga, the spectral width of a LIMO is inversely pro- 
portional to the oscillation power, similar to the famous Schawlow-Townes formula [23,35] for describing 
the spectral width, A viaser, of a laser: 


A l^laser — 


_ Ps 
27T ^~laser^ aser 


(28b) 


where p s = kv is the spontaneous emission noise density of the laser, Pi aS er is the laser oscillation power, 
and Ti aser is the decay time of the laser cavity. However, because, as will be shown in Section III.E, 
both P osc and px are functions of the photocurrent, the statement that the spectral width of a LIMO is 
inversely proportional to the oscillation power is valid only when thermal noise dominates in the oscillator 
at low photocurrent levels. 

From Eq. (28a), the quality factor, Q, of the oscillator is 


Q = 


fo. 


A Ifwhm 



where Qo is the quality factor of the loop delay line and is defined as 


(29) 


Qd = 2 t rf osc T 


(30) 


From Eq. (27), we easily obtain 
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(31a) 


SrfW) 


4 r 2 
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i/'l « 


AIfwhm 

2 


SflF(/0 


6 

( 27 t ) 2 ( t /') 2 


l/'l » 


AJfwhm 

2 


(31b) 


It can be shown [24] that for an oscillator with a phase fluctuation much less than unity, its power spectral 
density is equal to the sum of the single side-band phase noise density and the single side-band amplitude 
noise density. In most cases in which the amplitude fluctuation is much less than the phase fluctuation, 
the power spectral density is just the single side-band phase noise. Therefore, it is evident from Eq. (31b) 
that the phase noise of the LIMO decreases quadratically with the frequency offset, /'. For a fixed 
the phase noise decreases quadratically with the loop delay time. The larger the r, the smaller the phase 
noise. However, the phase noise cannot decrease to zero no matter how large r is, because at large enough 
r, Eqs. (27) and (31b) are not valid anymore. From Eq. (27), the minimum phase noise is Sjjj? ~ 6/4 at 
/' = l/2r. For the frequency offset, /', outside of the passband of the loop filter (where F(f f ) = 0), the 
phase noise is simply the noise-to-signal ratio, <5, as can be seen from Eq. (22). 

Equations (25), (27), and (31b) also indicate that the oscillator’s phase noise is independent of the 
oscillation frequency, f osc . This result is significant because it allows the generation of high-frequency and 
low phase- noise signals with the LIMO. The phase noise of a signal generated using frequency-multiplying 
methods generally increases quadratically with the frequency. 

E. The Noise-to-Signal Ratio 

As mentioned before, the total noise density input to the oscillator is the sum of the thermal noise, 
P thermal = 4ksT{NF)\ the shot noise, p s hot — 2 eI p h,R\ and the laser’s relative intensity noise (RIN), 
Prin = N RIN I* h R, densities [25,26]: 


Pn = ikBT(NF) + 2 eI p hR + N R iNl^ h R 


(32) 


where ks is Boltzmann’s constant, T is the ambient temperature, NF is the noise factor of the RE 
amplifier, e is the electron charge, I p h is the photocurrent across the load resistor of the photodetector, 
and Nrin is the RIN of the pump laser. 

From Eqs. (26) and (32), one can see that if the thermal noise is dominant, then 6 is inversely 
proportional to the oscillating power, Pose of the oscillator. In general, P osc is a function of photocurrent, 
Iph, and amplifier gain, G A , as determined by Eqs. (18a), (18b), and (18c), and the noise-to-signal ratio 
from Eq. (26) is thus, 


IG^I 2 4kT(NF) + 2eI ph R 4- N RIN I* h R 
“ 1 - 1/|G S | 4r ? 2 co S 2(7r Vb/V^D^R [ 

111 deriving Eq. (33), Eqs. (4) and (18b) are used. From Eq. (33), one can see that 6 is a nonlinear 
function of the small signal gain of the oscillator. As shown in Fig. 5(a), it reaches the minimum value 
at | G s { — 3/2: 


^min — 


4kT(NF) + 2 eI ph R + N RIN P ph R 
(16/27 )P ph R 


(34) 
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Fig. 5. Calculations of the LIMO's input noise-to-signal ratio (a) as 
a function of small signal gain, IG S I, showing a minimum value at 
lGgl= 1.5 and (b) as a function of photocurrent for different values of 
the laser s RIN noise. 


where rj = 1 and cos^Vjg/V^) = 1 are assumed. The oscillation amplitude at IG 5 I = 3/2 can be obtained 
from Eq. (18b) as 


_ 2\Z2Vtc _ p r 0 T/ 

Kosc — ^ 0.52 Vv 

(?rV3) 


and the corresponding BF power is 


P — 

J O.SC 


4v; 2 _ iop^ dB 

(3tt 2 P) “ 3 


(35a) 


(35b) 


where P^ dB is the input 1 -dB compression power of the E/O modulator [26,27]. From Eq. (35b), one can 
conclude that, in order to have minimum noise, the oscillation power measured at the input of the E/O 
modulator should be 5 dB above the 1 -dB compression power of the modulator. Equation (35a) indicates 
that the noise of the oscillator is at minimum when the oscillating amplitude is roughly half of V n or the 
voltage in the oscillator is varying between the peak and the trough of the sinusoidal transmission curve 
of the E/O modulator. This makes sense because the modulator has its minimum sensitivity to voltage 
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variations at the maximum and the minimum of the transmission curve, and the most likely cause of 
voltage variations in a LIMO is the noise in the loop. 

It is evident from Eq. (34) that the higher the photocurrent, the less the noise-to-signal ratio of the 
oscillator until it flattens out at the laser’s RIN level. Therefore, the ultimate noise-to-signal ratio of 
a LIMO is limited by the pump laser’s RIN. If the RIN of the pump laser in a LIMO is -160 dB/Hz, 
the ultimate noise-to-signal ratio of the oscillator is also —160 dB/Hz, and the signal-to-noise ratio is 
160 dB/Hz. Figure 5(b) shows the noise-to-signal ratio, <5, as a function of photocurrent, / P / M for different 
RIN levels. In the plot, the small signal gain, Gs , is chosen to be a constant of 1.5, which implies that 
when I p f h is increased, the amplifier gain, G A , must be decreased to keep Gs constant. From the figure, 
one can easily see that 8 decreases quadratically with I p h at small I p h and flattens out at the RIN level 
at large I ph . 


F. Effects of Amplifier Nonlinearity 

In the discussions above, we have assumed that the nonlinear distortion of a signal from the E/O 
modulator is more severe than from the amplifier (if any) used in the oscillator, so that the oscillation 
amplitude or power is limited by the nonlinear response of the E/O modulator. Using an engineering 
term, this simply means that the output 1-dB compression power of the amplifier is much larger than the 
input 1-dB compression power of the E/O modulator [26]. 

For cases in which the output 1-dB compression power of the amplifier is less than the input 1-dB 
compression power of the modulator, the nonlinearity of the amplifier will limit the oscillation amplitude, 
V osc (or power P osc ), of the oscillator, resulting in an oscillation amplitude less than that given by 
Eqs. (18a), (18b), and (18c). The exact relation between the oscillation amplitude (or power) and the 
small signal gain, Gs, can be determined using the same linearization procedure as that used for obtaining 
Eqs. (18a), (18b), and (18c) if the nonlinear response function of the amplifier is known. However, all 
the equations in Section III.D for describing the spectrum of the oscillator are still valid, provided that 
the oscillation power in those equations is determined by the nonlinearity of the amplifier. For a high 
enough small signal gain, Gs, the oscillation power is approximately a few dB above the output 1-dB 
compression power of the amplifier. 

In all the experiments below, the output 1-dB compression power of the amplifiers chosen is much 
larger than the input 1-dB compression power of the modulator, so that the oscillation power is limited 
by the modulator. 

IV. Experiments 

A. Amplitude Versus Open-Loop Gain 

We performed measurements to test the level of agreement of the theory described above with exper- 
imental results. In all of our experiments, we used a highly stable diode-pumped Nd:YAG ring laser [35] 
with a built-in RIN reduction circuit [36] to pump the LIMO. The experimental setup for measuring the 
oscillation amplitude as a function of the open-loop gain is shown in Fig. 6(a). Here an RF switch was 
used to open and close the loop. While the loop was open, an RF signal from a signal generator with the 
same frequency as the oscillator was injected into the E/O modulator. The amplitudes of the injected 
signal and the output signal from the loop were measured with an oscilloscope to obtain the open-loop 
gam, which was the ratio of the output amplitude to the injected signal amplitude. The open-loop gain 
was varied by changing the bias voltage of the E/O modulator, by attenuating the optical power of the 
loop, or by using a. variable RF attenuator after the photodetector, as indicated by Eq. (4). When closing 
the loop, the amplitude of the oscillation was conveniently measured using the same oscilloscope. We 
measured the oscillation amplitudes of the LIMO for different open-loop gains at an oscillation frequency 
of 100 MHz, and the data obtained are plotted in Fig. 4(b). It is evident that the experimental data 
agreed well with our theoretical predictions. 
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Fig. 6. Experimental setups for measuring (a) the oscillation amplitude of a LIMO as a function of the 
small signal gain and (b) the phase noise of a LIMO using the frequency discrimination method. 


B. Phase Noise Measurement Setup 

We used a frequency discriminator method [28] to measure the phase noise of the LIMO, and the 
experimental setup is shown in Fig. 6(b). The advantage of this method is that it does not require a 
frequency reference and, hence, can be used to measure an oscillator of any frequency. Using a microwave 
mixer in the experiment, we compared the phase of a signal from the electrical output port of the LIMO 
with its delayed replica from the optical output port. The length of the delay line is important because, 
the longer the delay line, the lower the frequency offset at which the phase noise can be accurately 
measured. On the other hand, if the delay line is too long, the accuracy of the phase noise at a higher 
frequency offset will suffer. The length of delay used in our experiment is 1 km, or 5 /is. Because of this 
delay, any frequency fluctuation of the LIMO will cause a voltage fluctuation at the output of the mixer. 
We measured the spectrum of this voltage fluctuation with a high dynamic-range spectrum analyzer and 
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transferred the spectral data to a computer. Finally, we converted this information into the phase noise 
spectrum of the LIMO according to the procedures given in [28]. In these experiments, the noise figure 
of the RF amplifier was 7 dB. 

C. Phase Noise as a Function of Offset Frequency and Loop Delay 

Figure 7(a) is the log versus log scale plot of the measured phase noise as a function of the frequency 
offset, /'. Each curve corresponds to a different loop delay time. The corresponding loop delays for curves 
1-5 are listed adjacent to each curve, and the corresponding oscillation powers are 16.33, 16, 15.67, 15.67, 
and 13.33 dBm, respectively. Curve fitting yields the following phase noise relations as a function of 
frequency offset /': -28.7 - 201og(/'), -34.84 - 201og(/'), -38.14 — 201og(/'), -40.61 - 201og(/'), and 
-50.45 - 20 log(/'). Clearly, the phase noise has a 20-dB-per-decade dependence with the frequency 
offset, in excellent agreement with the theoretical prediction of Eq. (31b). 

Figure 7(b) is the measured phase noise at 30 kHz from the center frequency as a function of loop delay 
time, with data points extracted from curves 1-5 of Fig. 7(a) and corrected to account for oscillation power 
differences. Because the loop delay is increased by adding more fiber segments, the open-loop gains of the 
oscillator with longer loops decrease as more segments are connected, causing the corresponding oscilla- 
tion power to decrease. From the results of Fig. 8, discussed below, the phase noise of the LIMO decreases 




Fig. 7. Single side-band phase noise of a LIMO measured at 800 MHz: (a) phase noise 
spectra at different loop delays and their fits to Eq. (31 b) and (b) phase noise at a 30-kHz 
offset from the center frequency as a function of loop delay time. 
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linearly with the oscillation power. To extrapolate the dependence of the phase noise on the loop delay 
only from Fig. 7(a), each data point in Fig. 7(b) is calibrated using the linear dependence of Fig. 8, while 
keeping the oscillation power for all data points at 16.33 mW. Again, the experimental data agree well 
with the theoretical prediction. 

D. Phase Noise as a Function of Oscillation Power 

We have also measured the phase noise spectrum of the LIMO as a function of oscillation power, with 
the results shown in Fig. 8. In that experiment, the loop delay of the LIMO was 0.06 /is, the noise figure 
of the RF amplifier was 7 dB, and the oscillation power was varied by changing the photocurrent, I v h, 
according to Eqs. (3), (4), (18a), (18b), and (18c). With this amplifier and the photocurrent level (1.8 mA 
- 2.7 mA), the thermal noise in the oscillator dominates. Recall that, in Eqs. (25), (26), and (27), the 
phase noise of a LIMO is shown to be inversely proportional to the oscillation power. This is true if the 
gain of the amplifier is kept constant and the photocurrent is low enough to ensure that the thermal noise 
is the dominant noise term. In Fig. 8(a), each curve is the measurement data of the phase noise spectrum 
corresponding to an oscillating power, and the curves in Fig. 8(b) are the fits of the data to Eq. (31b). 
Figure 8(c) is the phase noise of the LIMO at 10 kHz as a function of the oscillation power, extracted 
from the data of Fig. 8(b). The resulting linear dependence agrees well with the theoretical prediction of 
Eq. (31b). 

E. Phase Noise Independence of Oscillation Frequency 

To confirm our prediction that the phase noise of the LIMO is independent of the oscillation frequency, 
we measured the phase noise spectrum as a function of the oscillation frequency, and the result is shown 
in Fig. 9(a). In the experiment, we kept the loop length at 0.28 /zs and varied the oscillation frequency 
by changing the RF filter in the loop. The frequency was fine tuned using an RF line stretcher. It is 
evident from Fig. 9(a) that all phase noise curves at frequencies 100, 300, 700, and 800 MHz overlap with 
one another, indicating a good agreement with the theory. Figure 9(b) is a plot of the phase noise data 
at 10 kHz as a function of the frequency. As predicted, it is a flat line, in contrast with the case when 
a frequency multiplier is used to obtain higher frequencies. This result is significant because it confirms 
that the LIMO can be used to generate high-frequency signals up to 75 GHz with a much lower phase 
noise than that which can be attained with frequency-multiplying techniques. 


V. Summary 

We have introduced a high-frequency, high-stability, high spectral purity, widely tunable electro-optic 
oscillator, which we have termed a LIMO. The high stability and spectral purity of the LIMO result from 
the extremely low energy storage loss realization obtained with a long optical fiber. The optical fiber 
is also virtually free of any frequency-dependent loss, resulting in the same long storage time and high 
spectral purity signals for both low- and high-frequency oscillation. On the other hand, the oscillation 
frequency of the LIMO is limited only by the speed of the modulator, which at the present can be as 
high as 75 GHz. As yet another unique feature, the output of the LIMO may be obtained directly as 
microwave signals or as intensity modulations on an optical carrier for easy interface with optical systems. 

We also analyzed the performance of this oscillator by deriving expressions for the oscillation threshold, 
Eq. (5); oscillation amplitude, Eqs. (18a), (18b), and (18c); and oscillation frequency, Eqs. (19a) and (19b). 
These results agree quite well with experimental data obtained with laboratory versions of the LIMO. 

We also derived the expression for the spectrum of the output of the LIMO and showed that it has a 
Lorenzian line shape, given by Eqs. (25) and (27). The spectral width of the output signal was found to 
be inversely proportional to the square of the loop delay time, given by Eqs. (28a) and (28b). In addition, 
at low optical pumping levels where thermal noise dominates, the spectral line width was found to be 
inversely proportional to the oscillation power of the oscillator, similar to the Schawlow -Townes formula 
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Fig. 9. Single side-band phase noise measurements of the LIMO at different oscillation 
frequencies: (a) phase noise spectra and (b) phase noise at a 10-kHz offset frequency as 
a function of osciiiation frequency, extracted from (a). The loop delay for the 
measurements is 0.28 \xs. 


describing the spectral width of a laser. Since increasing the optical pump power increases the oscillation 
power, in this regime the line width of the LIMO decreases as the optical pump power increases. On the 
other hand, at high pump powers where the pump laser's relative intensity noise dominates, the spectral 
width approaches a minimum value determined by the laser’s RIN noise, as given by Eqs. (25), (27), and 
(34). We measured the phase noise spectrum of the LIMO and verified our theoretical findings. 


It is important to note that the analysis performed here was for the specific case of the LIMO with 
a Mach Zehnder electro-optic modulator. Other modulation schemes, such as with electro-absorptive 
modulators or by direct modulation of semiconductor lasers [29], will also lead to signals with charac- 
teristics similar to those obtained in this work. For these examples, the theoretical approach developed 
above is still applicable after suitable modifications. The major change required in the analysis is the 
replacement of Eq. (1), which describes the transmission characteristics of a Mach-Zehnder modulator, 
with the appropriate equation for the specific modulation scheme. All other equations can then be derived 
in the same way as described in the theory. 
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Because of its unique properties, the LIMO may be used in a number of applications. As a voltage- 
controlled oscillator (VCO) [12], it can perform all the VCO functions in both electronic and photonic 
applications, including generating, tracking, and cleaning RF carriers. The LIMO has the unique property 
of actually amplifying injected signals [12] and, thus, may be used in high-frequency carrier regeneration 
and signal amplification. Other important potential applications of the LIMO include high-speed clock 
recovery [30,31], comb and pulse generation [12], high-gain frequency multiplication, and photonic signal 
upconversion and downconversion [32] in photonic RF systems [33,34]. 
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The performance of a scheme proposed for automated routine monitoring of 
deep-space missions is presented. The scheme uses four different tones (sinusoids) 
transmitted from the spacecraft (S/C) to a ground station with the positive iden- 
tification of each of them used to indicate different states of the S/C. Performance 
is measured in terms of detection probability versus false alarm probability with 
detection signal-to-noise ratio as a parameter. The cases where the phase of the 
received tone is unknown and where both the phase and frequency of the received 
tone are unknown are treated separately. The decision rules proposed for detect- 
ing the tones are formulated from average-likelihood ratio and maximum-likelihood 
ratio tests, the former resulting in optimum receiver structures. 


I. Introduction 

It has been proposed that automated routine monitoring of deep-space missions be provided by trans- 
mitting one out of n (typically n = 4) different subcarriers (tones) from the spacecraft (S/C) and then 
using a small automated terminal (for example, a 6-m low Earth orbiter terminal (LEO-T)-class) ground 
station to detect the presence or absence of each possible tone. The positive identification of each of 
the tones at the receiver will indicate different stages of the S/C, for example, S/C healthy, S/C needs 
help, S/C is going to transmit telemetry, etc. Since each of these tones is transmitted from the S/C to 
the ground over an additive white Gaussian noise (AWGN) channel along with the added possibility of 
Doppler distortion, the above-mentioned detection problem to be solved at the receiver can be formulated 
as a binary hypotheses test of signal plus noise versus noise only. In the most general ease, the signal 
is modeled as a constant power sinusoid with unknown [i.e., uniformly distributed on (— 7r,7r)] phase 
and unknown (i.e., uniformly distributed in some interval {fi,//) governed by the amount of Doppler) 
frequency. 

The optimum solution to problems of this nature is based upon maximum-likelihood (ML) consid- 
erations of the type discussed in VanTrees [1]. In particular, the solution takes the form of a binary 
hypothesis likelihood ratio test against a threshold whose value depends on the specified false alarm and 
detection probabilities, the available signal power- to- noise spectral density ratio, and the duration of 
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the observation of the hypotheses. We shall see that there are, in principle, two detection/estimation 
philosophies suggested by the ML approach, corresponding respectively to what is commonly known as 
noncoherent detection, wherein no attempt is made to estimate the unknown parameters prior to detec- 
tion, and pseudocoherent detection, wherein an attempt is made to first estimate the parameters (using 
an ML approach) and then to use these estimates to aid in the detection process [2]. Since there appears 
to be some question about which is the better approach, we shall consider both approaches, discuss their 
philosophical differences, and compare their performances. 

This article is organized in two parts. In Part 1, we consider the problem of optimally detecting 
a sinusoidal signal of known amplitude (power) and frequency but of unknown phase [i.e. , uniformly 
distributed on (-7r,7r)] transmitted from a S/C to the ground over an AWGN channel. In so far as 
the optimum receiver design is concerned, the problem will be formulated as a binary hypothesis test of 
signal plus noise versus noise only with a single unknown parameter (i.e., carrier phase). In Part 2, we 
consider the added possibility of Doppler distortion, which produces an uncertainty in the received carrier 
frequency. Once again, the problem can be formulated as a binary hypothesis test of signal plus noise 
versus noise only, where now the signal is modeled as a constant power sinusoid with unknown phase and 
unknown frequency. Unfortunately, however, the theory for this case is not a s well developed in [1] as for 
the case where frequency is known. Nevertheless, other researchers [3-6] have examined problems of this 
type in the context of frequency-hopped (FH) or direct sequence (DS) spread spectrum communication 
systems, and we shall make use of their results wherever appropriate. 


Part 1. Known Frequency and Unknown Phase 
II. The Average-Likelihood Ratio Test 

A. Derivation of the Optimum Decision Rule and Associated Receiver Structure 

Consider the transmission of a fixed (known) amplitude sinusoid with known frequency and unknown 
phase over an AWGN channel. As such, the received signal can be modeled over the interval of observation 
0 < t < T as corresponding to either of two hypotheses, namely, 


r(t) — s(t, 9) + n(t) = V2P cos(uj c t + 0) + n(t) 


(la) 


when indeed the signal was sent (hypothesis H 1 ) or 


r(t) = n(t) 


(lb) 


when the signal was not sent (hypothesis H 0 ). In Eq. (la), P,u) c respectively denote the known signal 
power and radian carrier frequency, and 6 denotes the unknown carrier phase assumed to be uniformly 
distributed in the interval (— 7 r, 7 r). Also, n(t) denotes the AWGN with single-sided power spectral density 
No W/Hz. 

The optimum detection of a signal transmitted over an AWGN channel is the solution to the problem 
of finding the likelihood ratio (LR), defined as the ratio of the conditional probability density functions 
(pdf’s) of the received signal under the two hypotheses, namely, 


A (r(t)) - 


P(r(QlJTi) 

p(r(t)\H 0 ) 


( 2 ) 


and then comparing this ratio to a suitable chosen threshold to decide between H 1 and Ho , i.e., 
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( 3 ) 


Hi 

A (r(t)) > V 

H o 


In the case where all parameters of the signal are known, the evaluation of the numerator and denominator 
of Eq. (2) is straightforward, namely, 


p(r(£)|tfi) 


exp 


_1_ 

Wo 


T 


J (r(t) - s(t)) 2 dt 
0 


p(r(*)|tf 0 ) 



( 4 ) 


When the signal has an unknown parameter, e.g., the phase 0 , then to compute the numerator of Eq. (3), 
we must first condition the pdf p(r(t)\Hi) on the unknown parameter (9) and then average over this 
parameter, i.e., 


P 


(r(t) \Hi) = f 

J —IT 


p(r{t)\Hi,6)p e {6)d0 


( 5 ) 


where po(0) denotes the pdf of the unknown parameter 0. In our situation, the phase is assumed to be 
completely unknown and, thus, po(0) is a uniform distribution. Also note that this conditioning on the 
unknown parameter is now necessary in the denominator of Eq. (3) since the signal does not explicitly 
appear in p(r(t)\Ho) [see Eq. (4)]. Hence, combining Eqs. (3) through (5), the average-likelihood ratio 
(ALR) 2 becomes 


A(r(*)) - 


^ P(r(tm , D)dO S S’, “P {- ^ So 'MO - s(t, de 

p(r(tVH„) 1 f 1 r r 


J 

r pt) 

{ 1 . 

r j 

f 2 f T ( 

= exp <| 

{ ^oJ 

[ 2tt J 

1 exp < 

{No Jo r 

J 

r pt] 

{ 1 

r j 

{ 2 V 2 F r 

= exp j 

{ N 0 ] 

1 27 rj 

' exp < 

0 

£ 

✓ 


VtsN o 
r(t)s(t,9)dt ^ d6 


i 


(6) 


2 We shall refer to this formulation as an average-likelihood ratio (ALR) test to distinguish it from another (in general, less 
optimum) approach to be discussed shortly, in which a best (maximum-likelihood) estimate of the unknown parameter is 
obtained first and then used in the detection process. We shall refer to the latter approach as a maximum - likelihood ratio 
(MLR) test. This vernacular is not standard in the literature. What is important to understand here is that the words 
average and maximum refer to the manner in which the unknown parameter is handled, i.e., the estimation part of the 
problem and not the manner in which the decision on the hypothesis is made, i.e., the detection part of the problem. We 
shall be more explicit and mathematically precise about these differences later on in the article. 
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In arriving at the final result in Eq. (6), we have noted that the term exp j — / Q r r 2 (t)dt/No j is common 
to both the numerator and denominator and, thus, cancels, and also that 


exp \-'ki s2mdt }" exp {~^i c«v t < + s)4=»p{-f} 


( 7 ) 


assuming uj c T >> 1, as is typically the case. Defining the in-phase (I) and quadrature (Q) correlations 


[ r(t)y / 2 
Jo 


cos uj c tdt 


f r(t)V 2 
Jo 


sin uirtdt 


then Eq. (6) can be rewritten as 


A W f» = «P } £ /__ «P { ^ i “»(« + M } <» = <*P f - ^ } /o t 


(8) 


where 


l ± v'lfTif 


A _ 1 L s 
a = tan — 

-C«c 


( 9 ) 


Comparing A(r(t)) to a threshold rj is equivalent to comparing In A(r(t)) to In?;. Thus, taking the natural 
logarithm of Eq. (8), we obtain the equivalent decision rule 


Ho 


( 10 ) 


Finally, since In J 0 (x) is a monotonic function of its argument, x, and since PT/Nq can be absorbed into 
the decision threshold, then the decision rule of Eq. (10) can be further simplified to 


H x 

^1 L > < 
N 0 L < * 

Ho 


( 11 ) 


or, equivalently, 
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( 12 ) 


H i 

r2 > ,2^0 A 

< ? 4P 7 

H 0 


i.e., the optimum decision of signal present versus signal absent is determined from a comparison of 
the output of a square-law envelope detector with a normalized threshold, 7, whose value is determined 
from the specifications on false alarm probability and detection probability (see the next section). An 
implementation of the decision rule in Eq. (12) is illustrated in Fig. 1. 



Fig. 1. Average-likelihood (noncoherent) detector for detection of a single sinusoidal 
tone with known frequency and unknown phase in AWGN. 


B. Performance (Receiver Operating Characteristic) 

The performance of the receiver in Fig. 1 is described in terms of its false alarm probability ( Pf ), 
defined as the probability of deciding Hi (signal is present) when indeed Ho is true (signal is absent), and 
its probability of detection (Pd), defined as the probability of deciding Ho (signal is absent) when indeed 
Hi is true (signal is present). These probabilities are readily computed from knowledge of the first and 
second moments of the Gaussian random variables L c and L s [see Eq. (8)] under the two hypotheses, 
namely, 


Ho : E{L C } = E{L S } = 0 


var {L c } = var {L s } 


N 0 T 


H x : E{L c \e} = JPT cos 9 

E{L s \e } = -VPT sin 6> 


var {L c } = var {L s } = 


JVnT 


(13) 


To compute P F , we observe that, under hypotheses H 0 , L is a Rayleigh random variable (L 2 is a central 
chi-squared random variable). Thus, 
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* - pr[H ' w - priL 2 > - r /" (- 4 =) ^ 

= /j^ 2fiexp( - R2, " = “ p (-^) <i4) 

Similarly, we observe that, under hypothesis # 1 , L is a Rician random variable (L 2 is a noncentral 
chi-squared random variable). Thus, 


P D - - 1M* > •**.> - f (-£^1) /. (^) it 


/? 2 = (£{L C |0}) J + (£{L S |0}) 2 = PT 2 


> (15) 


J \j2~r/N,iT \ Z 


I 0 (Rd)dR = Q ( d, 


27 

AT n T 


where 


2 a 2PT = 2 E 
No ~ No 

is the detection signal-to-noise ratio (SNR) and Q(a,) 3) is the Marcum Q-function defined by [1]: 

t°° ( z 1 + q 2 \ 

Q(a,P)~ I 2 exp f j I 0 (az)dz 


(16) 


(17) 


Combining Eqs. (14) and (15) and eliminating the normalized detection threshold, one obtains the receiver 
operating characteristic (ROC) given by 


P D = Q (d, y/—2 In Pp) (18) 

which is illustrated in Fig. 2 for several values of the parameter d (or, equivalently, E/No). Alternatively, 
Pd I s plotted versus d 2 with Pp as a parameter in Fig. 3. 


III. The Maximum-Likelihood Ratio Test 

A. Derivation of the Optimum Decision Rule and Associated Receiver Structure 

Although the exact evaluation of the numerator of the likelihood ratio in Eq. (2), i.c., p(r(t)\H\) is 
obtained from the law of conditional probability as described by Eq. (5), namely, conditioning on the 
unknown parameter and averaging its distribution, it is also possible to approximate this numerator 
by first finding the ML estimate of the unknown parameter and then substituting this value into the 
conditional probability p(r(t)\H\,9). That is, we approximate p(r(t)\Hi) by p ^r(£)|//i , * in which 

case the likelihood ratio test (now referred to as the maximum-likelihood ratio (MLR) test) becomes 
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Fig. 2. ROC: frequency known and phase unknown. 


A (r(t)) ~ 


H i 

p{r{t)\H u 9 ML ) > 
p(r{t)\H 0 ) < 


Ho 


(19) 


We refer to this approach of first optimally estimating the phase and then using this estimate to aid 
the detection process as pseudocoherent detection. It is important at this point to emphasize that in the 
general context of problems of this type, i.e., detection of signals with completely unknown parameters, 
the performance of a receiver derived from MLR considerations (e.g., a pseudocoherent receiver) is never 
better than the performance of the ALR receiver (e.g., a noncoherent receiver ), which is indeed optimum 
under the assumed conditions. Thus, at best , one could hope that the MLR receiver would perform equally 
well as does the ALR receiver. In the next section, we shall indeed reveal the extent to which this equality 
in performance can be achieved for the problem at hand. First, however, let us derive the ML estimate 
of phase, namely, 9 ml, to be used in approximating the numerator of the likelihood ratio. 


The ML estimate of 6 is defined as 


9 m r = max 


p{r(t)\Hi,6) 
V p(r(t)\H 0 ) 


( 20 ) 


Using Eq. (4) in Eq. (20), it is straightforward to show that 


0 ml = max 
e 


exp {^oi> ^m*’*)*! = 


r 2V2 p p 

max exp < — - — / r(t) cos (a ) c t + 0) 
Q [ J o 


dt 


= max exp < ^^-Lcos(0 -F a) 
e iV 0 


( 21 ) 
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SINGLE TONE 
DOUBLE TONE 




Fig. 3. Detection probability {Pq) versus detection SNR (d 2 ): 
(a) frequency known and phase unknown and (b) frequency 
known and phase unknown (expanded view). 


where the envelope, L, and the phase, a, are defined by Eq. (9) together with Eq. (8). Since L is 
positive and independent of 0, then the maximization required in Eq. (21) is achieved when the argument 
of tiie cosine function is equal to zero. Thus, 


L = 


( 22 ) 


An implementation of this ML estimator of the unknown channel phase is illustrated in Fig. 4. Also 
illustrated in Fig. 4 is the pseudocoherent detector that employs this ML phase estimator, which can be 
obtained by taking the natural logarithm of Eq. (19). We now find the decision rule based on the MLR 
test in Eq. (19) and compare it with that of the previously discussed ALR test. Using Eq. (22) in Eq. (19) 
gives, by analogy with Eq. (8), 




MAXIMUM-LIKELIHOOD PHASE ESTIMATOR 


Fig. 4. Maximum-likelihood phase estimator and pseudocoherent detector. 


r pt) f 2\/P /- ( PT) f 2\/P 

A(r(*)) exp exp Lcos ( 0 ml + «J | = exp / exp 


(23) 


Taking the natural logarithm of Eq. (23), we then have, by analogy with Eq. (10), 


/- Hl 

2 VP r > . PT 

— — L _ In 77 4- — 

Nq < A^o 

^0 


(24) 


Since, as previously noted, the term PT /No can be absorbed into the decision threshold, then an equiv- 
alent test to Eq. (24) is 


2 sfp 
No 


L 





Ho 


(25) 


which, is identical to Eq. (11)! Thus, we conclude that in this particular circumstance, the MLR test 
(pseudocoherent receiver ) and the ALR test (noncoherent receiver) are identical Hence, the performance 
of the pseudocoherent receiver is also described by Figs. 2 and 3. It is to be emphasized again that 
the equivalence found here between ALR and MLR receivers is not typical and applies only in this very 
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special case of the detection of a signal with known frequency and unknown phase. More often than 
not, the receiver derived from the MLR approach will have an inferior performance to the optimum one 
derived from the ALR approach. 


IV. A More Precise Formulation of the Problem 

In reality, the subcarriers that are transmitted to indicate the status of the S /C are continuous square 
waves that biphase modulate the carrier. Thus, denoting the carrier radian frequency and phase by u c 
and 9 C (previously called 9 ), respectively, and the square-wave subcarrier radian frequency and phase by 
u sc and 0 SC , respectively, then the received signal analogous to Eq. (la) is given by 

r(t) = s(t , 9 C) 9 SC ) + n{t) = \/2Psin (u c t + 9 C -b ^ Sq ( u sc t -b 0 SC )) + n(t) 


= a/ 2 P Sq (u > sc t + 0 SC ) COS (uf c t + 0c) + n W 


(26) 


Assuming that the harmonics with frequencies other than the sum and difference of u c and u sc are filtered 
out, then in so far as detection is concerned, we may consider the received signal to be 3 

r(t) = s (t,9 c ,9 sc )+n(t) = \[P {cos[(cl> c + u; sc )t + (0 C + 0 SC )\ + cos[(u> c - u sc )t + (6 C - 0 3 c)]}+™(<) (27) 


i.e., the problem is to detect the presence or absence of two tones in an AWGN background where both 
uj c and Use are assumed to be known but both 9 C and 9 SC are assumed to be completely unknown. For 
convenience of notation, we shall rewrite Eq. (27) as 


r(t) = s(£,0+,0_) + n(t) — y/P {cos[u+t 4- 0+] + cos [u-t + 0_]} -b n(t) 


(28) 


where 


U± — uJ c i U sc 

9± = 6 C ± 6 sc 


(29) 


At first glance, it might appear that, because the phases 9 C and 9 SC appear in the two signal tones as 
their sum and difference, the detection of these tones cannot be performed independently. Interestingly 
enough, 9+ = 9 C + 9 SC and 0_ = 9 C — 9 SC when reduced modulo 27T can be shown to be independent 
uniformly distributed random variables (see the Appendix). Thus, as we shall see shortly, the detection 
of two distinct sinusoidal tones with independent random phases in an AWGN background can be treated 
by a likelihood ratio approach analogous to that discussed in the previous section for a single tone in the 
same background. 


* In reality, the >/F amplitude factor in Eq. (27) should be (2\/2/7r)\/P = 0.9003\/P to account for the amplitude of the 
first harmonic in the square-wave subcarrier. For simplicity, we shall ignore this minor difference. 
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A. The ALR Test 

As discussed in Section II, the optimum decision rule is, in general, obtained by applying the average- 
likelihood approach, which in this case means averaging the conditional likelihood function over the two 
random phases 8 + and 0-. In particular, the conditional pdf of the received signal under hypothesis H\ 
is analogous to Eq. (5): 


P{r{t\H,) = /:/; p{r(t)\H u 8 + ,e-)po + ,e_{8 + ,8-)d8 + dd- = J J p(r{t)\H x , 8+, 8-)d8 + d8. 

(30) 


and, hence, the ALR becomes 


A(r{t)) 


2 ^) I-* J*„P{r{t)\H u d + ,8~)d8 + dd_ 


p(r(t)\H 0 ) 


= exp 


r pt] ( i \ 2 r (2 Vp r 

{ Nq J \27tJ /.-( No i 


( 2 Jp r T 1 

r(t) cos (uj-t 4- 9- )dt > dO. 


/ 7T f n / p /*‘J 

exp | ^ — J r(t) cos(uj+t + 0 + )d£ | 




Defining the I and Q correlations for the sum and difference frequencies by 


(31) 


Lc± 


l cos uj±tdt 


= [ r(t)y/2 c 

Jo 

rT 

L s ± = / r(£)\/2sinu;±td£ 


(32) 


then, the likelihood function of Eq. (30) can be rewritten as 




(33) 


where, analogous to Eq. (9), the envelopes corresponding to the upper and lower subcarrier tones are 
given by 


L ± “ + L 2 S± 


(34) 


Alternately, in terms of the log-likelihood function, we arrive at the decision rule 
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(35) 


In J 0 



+ In /o 



H i 


> 

< 


In rj + 


PT 

Ao 




Note that now, despite the fact that In Iq(x) is a monotonic function of its argument, x, we cannot directly 
simplify Eq. (35) to a form analogous to Eq. (11). Rather, to get such a form, one must approximate the 
ln/ 0 (x) function by its series and asymptotic forms for small and large arguments, namely, 


lnl 0 (x) = 


x 

— , small x 
4 

|x|, large x 


(36) 


Thus, for example, if we invoke the small argument approximation of the In Iq(x) function in Eq. (35), 
we get the decision rule (optimum for small SNR) 


Hi 

Ll + L\ 4 L 2 J 7 (37) 

Ho 

where 7 is again a normalized threshold [not necessarily equal to the one defined in Eq. (12)]. The 
decision rule in Eq. (37) suggests the ALR structure illustrated in Fig. 5, which is analogous to that given 
in Fig. 1. For the large argument approximation of the \nIo(x) function, the implementation of Fig. 5 
would require square root devices in each arm entering the final summer prior to the decision device. 

B. The MLR Test 

Let us now again compare the noncoherent two-tone detector derived from ALR considerations and 
specified by the decision rule of Eq. (35) to a pseudocoherent detector that can be derived from MLR 
considerations. In particular, consider the joint ML estimates of defined as 


@ m l+ 1 $ m L— = max 


p(r(t)\H u 6 + ,6_) 


*+.*- p(r(t)\H 0 ) 

which, because of the independence of 6 + and 0_, is determined as 


(38) 


$ M l+ ■> @ m L— = max 
9 + ,0 


f 2 \fp r T ) { 2 s/p f T 

x exp < ^ J r(t) cos (u-t + 6_)dt > exp < — - J r(t) cos(a j+t + 0 4 


)dt 


= »iax |^^L_cos(0_ +a_)|exp|^^L + cos(0 + +a+)| 
where a± are defined in terms of L±, analogous to Eq. (9). The solution to Eq. (39) is 


(39) 


®ML± = 


(40) 
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Fig. 5. Average-likelihood (noncoherent) detector for detection of a pair of independent 
sinusoidal tones with known frequencies and unknown phases in AWGN. 


which, when substituted in Eq. (39), gives 

A(r(())ae x p{-^}exp|^L + |exp|^L.|=exp{-^|exp|^£j (41) 

Taking the natural logarithm of Eq. (40), we get the decision rule 

^(£ + + L_)> M+f (42) 

Ho 

Comparing Eq. (42) with Eq. (35), we observe that, in the two-tone case, the MLR test (which would 
lead to a pseudocoherent form of detector analogous to Fig. 4) is not the same as the ALR test. However, 
using the large argument approximation of the lnio(x) function as given by Eq. (36), we see that the 
ALR and MLR tests once again become equivalent. In summary then, we observe that, for detection of 
a single tone in AWGN, the ALR (noncoherent) test and MLR (pseudocoherent) test are equivalent for 
all SNRs, whereas for the detection of a pair of equal power tones in AWGN, the ALR and MLR tests 
are equivalent only at sufficiently large SNR. 

C. Performance (Receiver Operating Characteristic) 

The performance of the low SNR receiver in Fig. 5 is, as before, described in terms of its false alarm 
probability ( Pp ) and its probability of detection (Pd)- These probabilities are readily computed from 
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knowledge of the first and second moments of the Gaussian random variables L c ± and L s ± [see Eq. (31)] 
under the two hypotheses, namely, 


Ho : E{L C ±} = E{L S ±} = 0 


var { L c ± } = var {L s± } = 


N 0 T 


H\ : E{L c± \0) = \ —T cos 0± 


E{L s± \e} = \/ -^T sin 9± 


var { L c ± } = var {T s ±} 


N 0 T 


(43) 


To compute P F , we observe as before that, under hypothesis Ho, L 2 is a central chi-squared random 
variable (now with two more degrees of freedom). Thus, 


P F = Pr {Hi\Hq} = Pr{L 2 > 7 ^ 0 } = T ^rexp(-r)dr = (l + exp 


(44) 


Similarly, we observe that, under hypothesis H\,L 2 is a noncentral chi-squared random variable (now 
with two more degrees of freedom). Thus, 


P D =Pv{H l \H i }=Pv{L 2 > 1 \H l }= I 


y/27/ NuT \ d 


R ( w ) exp 


R 2 + d 2 \ ^ ( Rd j dR = / 2 7 


N n T 


(45) 


where d 2 is the detection SNR defined as before [see Eq. (16)] and Q]^(a,j3) is the generalized Marcum 
Q-function defined by 


/ z \ ^ ~~ 1 - 4 - Ck^ \ 

Qm(<*,(3) = J z J exp y ) I M -\{otz)dz (46) 

Note that Qm{oi,P) can be obtained from Q(a,f3) = Q\(a,(3) by the relation [2, Appendix 5A] 

Q M (a,P) = Q(a,P) + ex p (°-±^ £ Ij(a0) 

Unfortunately, the normalized detection threshold cannot be explicitly eliminated in Eqs. (39) and (40) 
to give a closed-form expression for the receiver operating characteristic (ROC) analogous to Eq. (18). 
However, for any range of interest, the ROC can be determined numerically. Such numerical results are 
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superimposed on the single-tone detection in Figs. 2, 3(a), and 3(b). We observe that the performance 
penalty associated with using a pair of subcarrier tones each with half the total power relative to the 
full-power single carrier tone case is quite small, e.g., on the order of 0.4 dB or less for P F = 1CT 2 and 
on the order of 0.3 dB or less for P F = 10“ 4 . The degradation associated with the true optimum ALR 
scheme as described by the decision rule of Eq. (35) would be even smaller. Thus, the performance curves 
of the true optimum ALR scheme for two tones would lie between the solid and dashed curves in Figs. 2, 
3(a), and 3(b) since indeed these performance results cannot beat those corresponding to the single-tone 
case. Because of the small degradations involved, we choose not to simulate the true optimum case. 


Part 2. Unknown Frequency and Unknown Phase 
V. The Average-Likelihood Ratio Test 

A. Derivation of the Optimum Decision Rule and Associated Receiver Structure 

Consider the transmission of a fixed (known) amplitude sinusoid with unknown frequency and unknown 
phase over an AWGN channel. Analogous to Eq. (1), the received signal can be modeled over the interval 
of observation 0 < t < T as corresponding to either of two hypotheses, namely, 

7'(t) = s(t , 6) 4- n(t) = V2Pcos(ujt + 6) 4- n(t) (47a) 

when indeed the signal was sent (hypothesis Hi) or 


r(t) = n(t) 


(47b) 


when the signal wets not sent (hypothesis Hq). In addition to the previously defined parameters, in 
Eq. (47a), / = lu/2i r denotes the unknown carrier frequency assumed to be uniformly distributed in the 
interval (f c — B/2, f c + B / 2), where as before f c denotes the nominal carrier frequency (i.e., in the absence 
of Doppler). When the signal has two unknown parameters, e.g., the phase 0 and frequency /, then to 
compute the numerator of Eq. (3), we must first condition the pdf p(r(t)\H\) on both of the unknown 
parameters and then average over them, i.e., 

rfrA-B/2 rTT 

p(r(t) \Hi)=l / p(r(t)\H l) ej)p e (9)p f (f)dedf (48) 

Jfa-Bf 2 J-ir 

where Pe(&),Pf(f) respectively denote the pdf’s of the unknown parameters 6 and /. In our situation, 
the phase and frequency are assumed to be completely unknown, and thus pe(6) and p/(f) are uniform 
distributions. Hence, combining Eqs. (3) and (48), the average-likelihood ratio (ALR) becomes 


A(r(t)) 


1 

2nB 



p(r(t)\H u e,f)d6df 


p(r(t)\H 0 ) 


exp 


exp 


f_PT| J_ rf * +B ? 2 

I Wo / 27 tB Jf' -B/Z 

\ PT] 1 r^ +B / 2 

iNo/S J fc _ B/2 0 



4v/p 

No 


L(f)J df 



(49) 
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where 


L(f) 4 y/L\{f) + L\(f) 


with 


rT 

L c (/) = / r(t)\/2cos27Tftdt 

Jo 

L s (f) = J r(t)y/2 sin ft dt J 


(50) 


(51) 


It should be noted that L(/) is nothing more than the magnitude of the complex Fourier transform (FT) 
of r(t) in the interval 0 < t < T. If r(t) is band limited to W Hz, then for large WT , the real and 
imaginary components of this complex FT, namely, L c (f), L s {f) can be approximated by the discrete 
Fourier transforms (DFTs) 


- 2 WT 

LM = " ( 2 #) C “ 

n= 1 
2 WT 

i.(/) = Y, r ( 2 ^) siu ( 2 i /5w) 


(52) 


Comparing A {r(t)) to a threshold produces (after suitable normalization) the ALR decision rule 



(53) 


Since Eq. (53) is overly demanding to implement, one discretizes the frequency uncertainty interval into 
G = B/T~~ l - BT subintervals to each of which is associated a candidate frequency /*; i — 0, 1 , 2, * • • , G-\ 
located at its center. As such, the integration over the continuous uncertainty region in Eq. (53) is 
approximated by a discrete (Riemann) sum and, hence, the approximate decision rule becomes 


G - 1 


s> 



Hi 

> 

< ^ 
Ho 


(54) 


which has the implementation representation of Fig. 6. It is important to understand that spacing the 
frequencies 1 ~ 0, 1, 2, • • • , G — 1 by 1/T guarantees independence of the noise components that appear 
at the output of each spectral estimate channel. However, orthogonality of the signal components of these 
same outputs depends on the true value of the received frequency relative to the discretized frequencies 
f , \ i = 0, 1, 2, ■ • ■ , G— 1 assumed for implementation of the receiver. That is, if the true received frequency 
happens to fall 011 one of the fi s, then a signal component will appear only in the corresponding spectral 
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Fig. 6. Average-likelihood (noncoherent) detector for detection of a single sinusoidal tone with unknown frequency and 

unknown phase in AWGN. 


estimate channel, i.e., all other channels will be noise only. On the other hand, if the true received 
frequency falls somewhere between two of the ff s, then we have a loss of orthogonality in that a spillover 
of signal energy occurs in the neighboring spectral estimates. The worst-case spillover would occur when 
the true received frequency is midway between two of the ff s. 

We conclude this section by noting that a decision metric similar to Eq. (54) arises in the study of 
FH or DS/low probability of intercept (LPI) optimum ALR (noncoherent) detection [3-5], where in the 
FH case, /*; i — 0, 1, 2, ■ * • , <3 — 1 corresponds to the G possible frequencies that the transmitted signal 
can hop to and the detection is based on observation of a single hop of duration Th = T , and in the DS 
case G is the number of possible code sequences that can occur in the observation interval. Many of the 
results obtained from these works are directly applicable to the problem at hand. 

B. Performance 

It is tempting for large values of G (as is typically the case) to apply a central limit theorem argument 
to the left side of Eq. (11), i.e., approximate it as a Gaussian random variable in so far as computing the 
receiver operating characteristic associated with this decision rule [4]. Unfortunately, it was shown in [5] 
that following such an approach is very poor when compared with results obtained from simulation or 
numerical methods applied to the true decision rule of Eq. (11), even for values of G as large as 1000 or 
10,000. In fact, it is stated in [5] that G on the order of “ten thousands is not guaranteed to be large 
enough to validate the Gaussian approximation.” Thus, to obtain the true receiver performance, we too 
must resort to simulation and/or numerical methods, such as those suggested by Requicha [7], wherein 
the characteristic function and fast Fourier transforms (FFTs) are used to compute approximate values 
of the distribution function associated with the left-hand side of Eq. (11). More about this later on. 
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VI. The Maximum-Likelihood Ratio Test 

A. Derivation of the Optimum Decision Rule and Associated Receiver Structure 

Although the exact evaluation of the numerator of the likelihood ratio in Eq. (2), i.e., p(r(t)\H\) is 
obtained from the law of conditional probability as described by Eq. (5), namely, conditioning on the 
unknown parameters and averaging over their distribution, it is also possible to approximate this numer- 
ator by first finding the ML estimates of the unknown parameters and then substituting these values into 
the conditional probability p(r(t)\Hi, 0, /). That is, we approximate p(r(t)\Hi) by p(r(t)\Hi,Q ml, f ml), 
in which case the likelihood ratio test (now referred to as the ra^mmura-likelihood ratio (MLR) test) 
becomes 


A (r(t)) * 


P (r{t)\HuOb xl, /ml) 
p(r(t)\H 0 ) 


H i 
> 
< 
H 0 


where 


Oml, f ml = max 
vj 


p(r(t)\H 0 ) 


(55) 


(56) 


The maximization over 6 required in Eq. (56) can be performed identically to that in Section III [see 
Eq. (23)]: 


max 

9 


P(r(t)\H 

P(r(t)\H 0 ) 


= exp 




(57) 


where L(f) is as defined in Eq. (50). Thus, the optimum maximum a posteriori (MAP) decision rule 
becomes 


max exp 
/ 



H i 

> 

< 71 
H o 


(58) 


Since the exponential is a monotonic function of its argument, we have the equivalent decision rule 4 


H i 

max L(f) > y/j (59) 

Ho 

which results in a spectral m.axim,um form of receiver. Again, because of the excessive demand placed on 
the implementation by the need to evaluate Eq. (59) over a continuum of frequencies, we again quantize 
the frequency uncertainty interval into G = BT subintervals, each with an associated candidate frequency 


4 We define the normalized threshold equal to ^ to be consistent with the notation used in Part 1. In this way, when G 
is equated to unity, then our results obtained here will reduce to those given in Part 1 for the MLR decision rule. 
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G'Z. • 


fu i — 0, 1, 2, ■ ■ , G - 1 located at its center. Thus, the frequency continuous decision rule of Eq. (59) can 
be approximated by the decision rule 


H i 

maxL(/i) < V7 (60) 

Ho 

which suggests the receiver of Fig. 7. Here again, as with Eq. (54), the orthogonality of the spectral 
estimates is not guaranteed unless the frequency of the received signal falls on one of the /Vs. Also, since 
£(/;); * = 0, 1 , 2, • • • , G represents a uniform sampling of !(/), then in view of Eq. (52), we can implement 
Fig. 7 with FFT techniques. 



Fig. 7. Maximum-likelihood detector for detection of a single sinusoidal tone with unknown frequency and 

unknown phase in AWGN. 


B. Performance 

The performance of the MLR decision rule of Eq. (60) can be obtained analytically since the pdf of 
G independent random variables can be explicitly written in terms of the pdf’s of individual random 
variables, which in turn are obtained from the results in Part 1. The procedure is as follows. 

1. Best-Case Performance. Consider first the optimistic (best) case, where the actual received 
carrier frequency is indeed equal to one of the G frequencies, say fk used to approximately implement 
the optimum decision rule as per the discussion following Eq. (59). Under H\,G — 1 of the L(/i)’s are 
Rayleigh distributed with pdf [see Eq. (14)] 

= £>o ( 61 ) 
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and the single L(f l ) that is associated with the received signal carrier frequency, namely L(/fc), is Rician 
distributed with pdf [see Eq. (15)] 


PW „m = ^exp(-^)/„(^) ; L>0, s-pi* («> 

Let Pp denote the per frequency channel false alarm probability, i.e., 

which is independent of fi. Then, the overall false alarm probability, Pp, is given by 
P F = Pr | max L(/j) > V7l-^o} 

- 1- Pr {L(fo) < s/l-, L(fi) < y/i, • • • , L(fc-i) < y/l\Ho} 

G - 1 G - 1 

= 1 - n Pr W‘) ^ Vl\ H o} = 1 - n 0 - Pr W*) > yfr\H 0 }) 

1=0 1=0 

.l-(l-p?) G =l-(l-«p (64) 

Since, under Hi, G — 1 of the spectral estimates (i.e., the ones containing noise only) have the same pdf, 
namely Eq. (61), as under Hq , and one spectral estimate has the Rician pdf of Eq. (62), then the overall 
probability of detection, P d , is determined from 

G - 1 

P D = Pr | max L(f l ) > = 1- (1 - Pr{L(/;) > y/j\H],}) = 1-(1 - P F ) G 1 (1 - Pp) ( 65 ) 

1=0 

where Pb corresponds to the detection probability of the 
i.e., 


P' D = Pr {£,(/,)> s /j\H l } = Q^d. 
Substituting Eqs. (63) and (65) into Eq. (66) gives 

p D = (-^)y 

or, equivalently, the overall probability of miss, P F j, is 


single-frequency channel containing the signal, 


27 \ , 2 2PT 

NoTj ’ N 0 


( 66 ) 




( 67 ) 
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P " 4 l - Pc= (l-e x p(-A^)) 


(68) 



Note that, for G = 1, Eqs. (64) and (67) reduce, respectively, to Eqs. (14) and (15). 

The ROC can be determined by eliminating the normalized threshold between Eqs. (64) and (67), in 
which case one obtains 


Pd = 1 - (1 - Pf) {G ~ 1)/G (l-Q (d, y/-2]n(l-(l-P F )'/G)y) 


(69) 


2. Worst-Case Performance. The worst-case performance occurs when the actual received carrier 
frequency is indeed midway between two of the G frequencies used to approximately implement the 
optimum decision rule as per the discussion following Eq. (59). Under Hq, the false alarm performance is 
still described by Eq. (64). However, under H i, all G spectral estimates are now Rician distributed with 
pdf’s of the form in Eq. (62), namely, 


PLUi){L ) - 


2 

No T 


Lexp 



L > 0 


where the /Vs are determined as follows. Since [see Eq. (15)] 


(70) 


0i = (E {L e (fi)\9, f}f + (E {L a (fi)\e, f }f 


(71) 


then, assuming that the actual received carrier frequency, /, is situated midway between fa and fa+ 1 , 
which are separated by 1/T, i.e., / = f k -f 1/2T, Eq. (70) is evaluated as (for simplicity, we ignore the 
edge effects at the ends of the frequency uncertainty band) 



i — fc, k + 1 
i ^ k,k 4-1 


= pr 2 r t 

Finally then, analogous to Eq. (70), the detection probability would be given by 


(72) 


p D = i- n (i - Q ( r i 

i=0 ' ^ 



(73) 


which, in general, depends on /^, i.e., the location of / within the uncertainty band. 

It has been suggested in [5] that the two nearest spectral estimates (envelopes) to the frequency 
location of the received signal dominate the performance, i.e., the spillover effect of signal in the other 
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frequency slots can be ignored to a first-order approximation. When this is done, then, under H\, two of 
the spectral estimates are identically Rician distributed and the remaining G - 2 are identically Rayleigh 
distributed. In this case, Eq. (73) is replaced by an expression somewhat like Eq. (67), namely, 


Pd = 1 



(74) 


which is now independent of the frequency location of the signal. Combining Eqs. (64) and (74), the 
ROC is approximately given by 


P D = 1 - (1 - P F ) iG 


— 2 )/G 


1 - Q ((f) d * >/ _ 21n(l - (1 -Pf) 1/g 


(75) 


VII. A More Precise Formulation 

As discussed in Section IV of Part 1, the true transmitted signal corresponds to a sinusoidal carrier 
phase modulated by a square- wave subcarrier of radian frequency uj sc . At the receiver, the harmonics 
with frequencies other than the sum and difference of tu sc and lj c are filtered out, which means that in so 
far as detection is concerned, the received signal in the absence of frequency uncertainty can be modeled 
as 


r(t) = s(t,6 c , 6 sc ) + n(t) = y/P { cos [(cj c +c o sc )t -f (6 C + <9 SC )] + cos [(w c - co sc )t + ( 6 C - <9 SC )]} + n(t) 

(76) 


In the presence of frequency uncertainty due, for example, to Doppler shift, both the upper and lower 
frequency tones in Eq. (76) will be shifted from their nominal values with the higher-frequency tone 
experiencing a larger shift than that corresponding to the lower-frequency tone. If, however, the subcarrier 
frequency is much smaller than the carrier frequency, i.e., lj sc << u c , as is the case of interest, then for all 
practical purposes, one can associate the frequency uncertainty with the carrier as discussed in Section V. A 
and assume to a first-order approximation that both upper and lower frequency tones experience the same 
frequency shift. Stated another way, we can assume that, in so far as detection is concerned, we observe 
a pair of tones whose frequencies are unknown (but by the same amount), each in a band D Hz centered 
around its nominal value. Furthermore, the uncertainty band is assumed to be very narrow with respect 
to the subcarrier frequency, i.e., B « f sc . 

A. The ALR Test 

Analogous to what was done in Part 1, the conditional pdf of the received signal under hypothesis H\ 
is given by 


p(r(t) m 




p(r(t)\H u 6 + , 0_, / - / sc , / + f sc )d6+dd-df 


(77) 


whereupon the ALR becomes 
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A(,«» = ex P {^} 5 «t 


(78) 


In Eqs. (77) and (78), the spectral envelopes at the lower and upper tones are defined by 


L±(f) = 


( 79 ) 


together with 


L c ±{f) = [ r(t)V 2 cos[ 2 Tr(f ± f sc )t]dt 

Jo 

L s ±(f) = J r (t)V 2 sin [ 2 ir(f ± f sc )t]dt j 


Discretizing the integration interval results in the approximate decision rule 




X H ' 

L-(fi ) J 
' Ho 


(80) 


(81) 


where the spectral envelopes required in Eq. (81) are defined analogously to Eqs. (79) and (80), with the 
continuous random variable / replaced by the discrete random variable /*; i = 0, 1, • • • , G — 1. As was 
the case for the single-tone result in Section V.A, the performance (ROC) of the decision rule in Eq. (81) 
cannot be obtained analytically. 

B. The MLR Test 

Without going into great detail, it is straightforward to show (using the results of Section IV. B) that 
the MLR test analogous to Eq. (58) becomes 


T exp exp (V M/) ) exp (‘ W LM) ) < ” 

Ho 

or, equivalently, 

Hi 

ma x{L-(f) + L+(f)) > ^7 
Ho 

which has the discretized version 

Hi 

max(L_(/i) + £+(/)) < a/7 

l i 

Ho 


(82) 


(83) 


(84) 
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Unfortunately, the performance of the receiver that implements the decision rule of Eq. (84) also cannot 
be obtained analytically. 

VIII. Numerical Results 

Since the performance of none of the ALR optimum decision rules can be evaluated analytically and 
since the same is true for some of the MLR decision rules, a computer simulation of these metrics has 
been developed to numerically evaluate such performance. The results of such simulations are described 
as follows. Figure 8 is a sample illustration of the ROC for the case of a single tone with unknown phase 
and frequency (as described in Section V) and a detection SNR d 2 — 2PT/Nq = 6 dB. Both ALR and 
MLR cases are illustrated, corresponding, respectively, to the decision rules of Eqs. (54) and (60). Also, 
both the best- and worst-case input frequency scenarios are considered, corresponding, respectively, to the 
cases where the actual input frequency is indeed equal to one of the G frequencies used to approximately 
implement the decision rule and the case where the actual input frequency falls midway between any 
two of these G frequencies. Clearly, the actual system performance corresponding to an input frequency 
arbitrarily chosen in the uncertainty band will lie between these two performance bounds. We observe 
from the results in Fig. 8 that the difference between best- and worst-case performance is relatively small, 
as well as is the difference between the ALR (optimum) and MLR (suboptimum) decision rules. There 
is a significant difference, however, between the performance for G = 10 and G = 100, indicating the 
sensitivity of the performance degradation to a factor of 10 increase in frequency uncertainty. Also, 
comparing Fig. 8 with the analogous curve in Fig. 2, corresponding to the case of unknown phase but 
known frequency, we again see a rather significant degradation in performance when the frequency is 
unknown even by only a factor of 10 relative to the observation bandwidth (reciprocal of the observation 
time, T), i.e., G = 10. 



Fig. 8. ROC: frequency and phase unknown (single-tone) 
simulation results. 


As verification of the MLR simulation results, we present in Fig. 9 the analogous analytical results 
obtained from Eqs. (69) and (75). Recall that in arriving at Eq. (75) the assumption was made that 
t he energy spillover effect of the signal into the other frequency slots is dominated by the two adjacent 
ones. Thus, ignoring edge effects, it was not necessary to average over all possible worst-case (mid- 
way) input frequency positions. In the computer simulation, this assumption was not invoked, as the 
input frequency was allowed to occur midway between any two adjacent frequencies. Despite this analysis 
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approximation, however, comparison of the results in Figs. 8 and 9 reveals excellent agreement between 
analysis and simulation, i.e., the assumption of only adjacent signal energy spillover used to arrive at 
Eq. (75) has been justified. Also indicated in Fig. 8 is the analytical result corresponding to known phase 
and frequency (recall that this result is the same for both MLR and ALR) that allows a more direct 
assessment of the performance degradation due to lack of perfect frequency knowledge. 

Since the curves in Fig. 8 are drawn for a fixed value of detection SNR d 2 = 2PT/N 0) then assuming 
that P/N 0 is specified, this implies that the observation interval, T, is also held constant. Thus, changing 
the value of G = BT from 10 to 100 directly translates into a change by a factor of 10 in the frequency 
uncertainty region B , which accounts for the observed degradation in performance. Another interpretation 
of the numerical data can be obtained by again holding P/No fixed but observing the effect on system 
performance of increasing T for a fixed frequency uncertainty region B. This necessitates plotting the 
ROC with both d 2 and G increasing linearly with T. Such a plot for the ALR decision rule with best-case 
input frequency is illustrated in Fig. 10, where the ROC is plotted for values of G ~ 10, 20, 40, and 80 (T 
increasing by a factor of 2) and corresponding values d 2 = 6, 9, 12, 15 dB. To directly see the dependence 
of MLR system performance on detection SNR, Fig. 11 illustrates the behavior of detection probability, 
Pd, versus detection SNR, d 2 , for a fixed false alarm probability, Pp = 10 -2 , and values of G — 10 and 
100. These curves are obtained from numerical evaluation of the analytical results in Section VI. Since 
along any curve G is held fixed, one can interpret these results as keeping the frequency uncertainty band, 
R, and observation time, T, fixed and observing the change in performance as P/No is varied. 

The penalty associated with detecting a pair of subcarrier tones (each at half the total transmitted 
power) as opposed to a single carrier tone (at full transmitted power) is illustrated by the numerical 
results in Fig. 12. Here we plot the ROC for both the single- and double-tone cases for the ALR decision 
rule with best-case input frequency and a detection SNR equal to 6 dB. The results for the single-tone 
case are taken directly from Fig. 8. We observe a significant performance penalty associated with using a 
double- tone detection scheme. Figure 13 illustrates for the double- tone detection scheme results analogous 
to Fig. 10 for the single-tone detection scheme. Here again, by comparing the two figures, we observe a 
significant penalty associated with using a pair of equal half-power subcarrier tones rather than a single 
tone at full power. 



Fig. 10. ROC simulation results: frequency and phase unknown 
(single tone), ALR , best-case input frequency. 
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Fig. 13. ROC simulation results: frequency and phase unknown 
(double tone), ALR, best-case input frequency. 


Acknowledgment 

The authors would like to thank Mr. Van Snyder of the Supercomputing and 
Computational Mathematics Support Group for his advice on performing the nu- 
merical integrations required in many of the analytical results. 


References 

[1] H. L. VanTrees, Detection, Estimation, and Modulation Theory, Part 1 , New 
York: John Wiley & Sons, Inc., 1968. 

[2] M. K. Simon, S. M. Hinedi, and W. C. Lindsey, Digital Communication Tech- 
niques: Signal Design and Detection , Englewood Cliffs, New Jersey: Prentice- 
Hall, Inc., 1994. 

[3] A. Polydoros and C. L. Weber, “Optimal Detection Considerations for Low 
Probability of Intercept,” MI LOOM ’82 Conference Proceedings , Boston, Mas- 
sachusetts, pp. 2. 1-1-2. 1-5, October 1992. 

[4] A. Polydoros and K. T. Woo, “Wideband Spectral Detection of Unknown Fre- 
quency Signals,” presented at the International Symposium on Information The- 
ory, Brighton, England, June 1985. 

[5] C.-D. Chung and A. Polydoros, Multi-Hop FH/LPI Detection, Part II: Spectral 
Techniques , Technical Report CSI-87-09-01, University of Southern California, 
Los Angeles, California, September 1987. 

[6] A. Polydoros and C. L. Nikias, “Detection of Unknown-Frequency Sinusoids in 
Noise: Spectral Versus Correlation Detection Rules,” IEEE ASSP , vol. 35, no. G, 
pp. 897-900, June 1987. 

[7] A. A. G. Requicha, "Direct Computation of Distribution Function From Charac- 
teristic Functions Using the Fast Fourier Transform,” Proceedings of the IEEE , 
vol. 58, no. 7, pp. 1154-1155, July 1970. 


96 




Appendix 

On the Independence of the Sum of Difference of Two Uniformly 
Distributed Random Variables Modulo 2n 


Consider two independent random phases 9 A and 0 B that are each uniformly distributed in the semi- 
closed interval [— 7T, 7 t). Define the sum and difference of these two random variables by 


o' + =o A +e B ' 
o'_ =e A -e B , 


and the modulo 2n versions of these random variables by 


9+ - {9'+) mod 27r - (9a + 9b ), nod 2x 
«- = (*-)„* 2, “ («A - M mo d 2. 


(A-l) 


(A-2) 


The probability density functions (pdf’s) of 0' + and 6'_ are triangular in the semiclosed interval [— 27 t, 27 r), 
i.e., they are the convolutions of two uniform pdf’s, whereas the pdf’s of their modulo 27 t reduced versions, 
0+ and 0_, are once again uniformly distributed in [-7r,7r) (see Fig. A-l). We would now like to show 
that 0+ and 6 _ are indeed independent random variables. To do this, we shall show that the conditional 
pdf P0_(9-\O + ) satisfies p0_(0_|0+) = pe_(9 _), i.e., it is a uniform distribution in [— 7r,7r). Similarly, it 
can be shown that p0 + (O+\O _) — p# + (0 + ). 

Let 9+ be any positive value in its region of definition, i.e., 0 < 0+ < 7r. Then, 0 A and 0 B are related 
as follows: 


A _ f -0A 4- 0 + - 2i r, -7T < 9 a < -7T + 

B “ \ ~0 A + 0 + , -7T + 0 + < 9 A < 7T 

From Eq. (A-l), we find that 

_ j 29 A - 0+ 4* 2?r, -7t < 9 a < -7r + 0+ 

^ 29 a - 0+, -7T + 9 + < 9 a < 7i 


(A-3) 


(A-4) 


Thus, from Eq. (A-4) and the fact that 9 A is uniform in the interval [ — 7r, 7 t), the conditional pdf p e >_ (0'_ |0-j-) 
appears as in Fig. A-2(a). Reducing 0+ modulo 2tt produces the conditional pdf p0_ (0_ |0 + ) as illustrated 
in Fig. A-2(b), i.e., a uniform distribution in the interval [-7r,7r) Q.E.D. 
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In this article, we design new turbo codes that can achieve near- Shannon-limit 
performance. The design criterion for random interleavers is based on maximizing 
the effective free distance of the turbo code, i.e., the minimum output weight of 
codewords due to weight-2 input sequences. An upper bound on the effective free 
distance of a turbo code is derived. This upper bound can be achieved if the feedback 
connection of convolutional codes uses primitive polynomials. We review multiple 
turbo codes (parallel concatenation of q convolutional codes), which increase the 
so-called “interleaving gain ” as q and the interleaver size increase, and a suitable 
decoder structure derived from an approximation to the maximum a posteriori 
probability decision rule. We develop new rate 1/3, 2/3, 3/4, and 4/5 constituent 
codes to be used in the turbo encoder structure. These codes, for from 2 to 32 
states, are designed by using primitive polynomials. The resulting turbo codes have 
rates b/n, b=Z, 2, 3, 4, and n=2, 3, 4, 5, 6 and include random interleavers for better 
asymptotic performance. These codes are suitable for deep-space communications 
with low throughput and for near-Earth communications where high throughput is 
desirable. The performance of these codes is within 1 dB of the Shannon limit at a 
bit-error rate of 10~ 6 for throughputs from 1/15 up to 4 bits/s/Hz. 


I. Introduction 

Coding theorists have traditionally attacked the problem of designing good codes by developing codes 
with a lot of structure, which lends itself to feasible decoders, although coding theory suggests that 
codes chosen “at random” should perform well if their block sizes are large enough. The challenge 
to find practical decoders for “almost” random, large codes has not been seriously considered until 
recently. Perhaps the most exciting and potentially important development in coding theory in recent 
years has been the dramatic announcement of “turbo codes” by Berrou et al. in 1993 [7]. The announced 
performance of these codes was so good that the initial reaction of the coding establishment was deep 
skepticism, but recently researchers around the world have been able to reproduce those results [15,19,8]. 
The introduction of turbo codes has opened a whole new way of looking at the problem of constructing 
good codes [5] and decoding them with low complexity [7,2]. 

Turbo codes achieve near-Shannon-limit error correction performance with relatively simple component 
codes and large interleavers. A required E^/Nq of 0.7 dB was reported for a bit-error rate (BER) of 10~ 5 
for a rate 1/2 turbo code [7]. Multiple turbo codes (parallel concatenation of q > 2 convolutional codes) 
and a suitable decoder structure derived from an approximation to the maximum a posteriori (MAP) 
probability decision rule were reported in [9]. In [9], we explained for the first time the turbo decoding 
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scheme for multiple codes and its relation to the optimum bit decision rule, and we found rate 1/4 turbo 
codes whose performance is within 0.8 dB of Shannon’s limit at BER=10~ 5 . 

In this article, we (1) design the best component codes for turbo codes of various rates by maximizing 
the “effective free distance of the turbo code,” i.e., the minimum output weight of codewords due to 
weight-2 input sequences; (2) describe a suitable trellis termination rule for b/n codes; (3) design low 
throughput turbo codes for power-limited channels (deep-space communications); and (4) design high- 
throughput turbo trellis-coded modulation for bandwidth-limited channels (near-Earth communications). 


II. Parallel Concatenation of Convolutional Codes 

The codes considered in this article consist of the parallel concatenation of multiple (q > 2) con- 
volutional codes with random interleavers (permutations) at the input of each encoder. This extends 
the original results on turbo codes reported in [7], which considered turbo codes formed from just two 
constituent codes and an overall rate of 1/2. 

Figure 1 provides an example of parallel concatenation of three convolutional codes. The encoder 
contains three recursive binary convolutional encoders with mi, m2, and m3 memory cells, respectively. 
In general, the three component encoders may be different and may even have different rates. The first 
component encoder operates directly (or through 7Ti ) on the information bit sequence u = («i,- * * ,un) 
of length TV, producing the two output sequences xo and xi. The second component encoder operates 
on a reordered sequence of information bits, u 2 , produced by a permuter (interleaver), 7T 2 , of length TV, 
and outputs the sequence x 2 . Similarly, subsequent component encoders operate on a reordered sequence 
of information bits. The interleaver is a pseudorandom block scrambler defined by a permutation of TV 
elements without repetitions: A complete block is read into the the interleaver and read out in a specified 
(fixed) random order. The same interleaver is used repeatedly for all subsequent blocks. 

Figure 1 shows an example where a rate r = 1/n = 1/4 code is generated by three component codes 
with memory 7 U\ = m2 = ?n 3 = m = 2, producing the outputs x 0 = u, xj = u gi/go, x 2 = u 2 • tfi/tfcu 
and x 3 = u 3 ■ gi/go (here 71*1 is assumed to be an identity, i.e., no permutation), where the generator 
polynomials go and gi have octal representation (7 ) oc t a i and (5) oc t a /, respectively. Note that various code 
rates can be obtained by proper puncturing of xi, x 2 , x 3 , and even x 0 (for an example, see Section V). 

We use the encoder in Fig. 1 to generate an (n(N + m), TV) block code, where the m tail bits of code 2 
and code 3 are not transmitted. Since the component encoders are recursive, it is not sufficient to set 
the last m information bits to zero in order to drive the encoder to the all-zero state, i.e., to terminate 
the trellis. The termination (tail) sequence depends 011 the state of each component encoder after TV 
bits, which makes it impossible to terminate all component encoders with rn predetermined tail bits. 
This issue, which had not been resolved in the original turbo code implementation, can be dealt with 
by applying a simple method described in [8] that is valid for any number of component codes. A more 
complicated method is described in [18]. 

A design for constituent convolutional codes, which are not necessarily optimum convolutional codes, 
was originally reported in [5] for rate 1/n codes. In this article, we extend those results to rate b/n 
codes. It was suggested (without proof) in [2] that good random codes are obtained if g (l is a primitive 
polynomial. This suggestion, used in [5] to obtain “good” rate 1/2 constituent codes, will be used in this 
article to obtain “good” rate 1/3, 2/3, 3/4, and 4/5 constituent codes. By “good” codes we mean codes 
with a maximum effective free distance d e j , those codes that maximize the minimum output weight for 
weight-2 input sequences, as discussed in [9], [13], and [5] (because this weight tends to dominate the 
performance characteristics over the region of interest). 
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III. Design of Constituent Encoders 

As discussed in the previous section, maximizing the weight of output codewords corresponding to 
weight-2 data sequences gives the best BER performance for a moderate bit signal-to-noise ratio (SNR) 
as the random interleaver size N gets large. In this region, the dominant term in the expression for bit 
error probability of a turbo code with q constituent encoders is 






2r t(g^ + 2 


where d p 2 is the minimum parity-weight (weight due to parity checks only) of the codewords at the 
output of the jth constituent code due to weight-2 data sequences, and 0 is a constant independent of 
N. Define dj ^ = d p 2 + 2 as the minimum output weight including parity and information bits, if the j th 
constituent code transmits the information (systematic) bits. Usually one constituent code transmits the 
information bits (j ~ 1), and the information bits of others are punctured. Define d e f = + 2 38 

the effective free distance of the turbo code and 1 /N q ~ l as the “interleaver’s gain.” We have the following 
bound on d 2 for any constituent code. 

Theorem 1 . For any r = b/(b 4- 1) recursive systematic convolutional encoder with generator matrix 


hi(D) - 

ho(D) 

h*{D) 

h*{D) 


h b (D ) 

ho(D ) . 
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where Ibxb is a 6 x b identity matrix, deg[/ii(£»)] < m, hi(D) ± h 0 (D), i = 1,2, •••,6, and h 0 (D) is a 
primitive polynomial of degree m, the following upper bound holds: 


om— 1 

<*5<l— J+2 


Proof. In the state diagram of any recursive systematic convolutional encoder with generator matrix 
G , there exist at least two nonoverlapping loops corresponding to all-zero input sequences. If ho(D) is a 
primitive polynomial, there are two loops: one corresponding to zero-input, zero-output sequences with 
branch length one, and the other corresponding to zero-input but nonzero-output sequences with branch 
length 2 m — 1, which is the period of maximal length (ML) linear feedback shift registers (LFSRs) [14] 
with degree m. The parity codeword weight of this loop is 2 m_1 , due to the balance property [14] of ML 
sequences. This weight depends only on the degree of the primitive polynomial and is independent of 
h t {D), due to the invariance to initial conditions of ML LFSR sequences. In general, the output of the 
encoder is a linear function of its input and current state. So, for any output we may consider, provided 
it depends on at least one component of the state and it is not ho(D), the weight of a zero-input loop is 
2 m “ 1 , by the shift-and-add property of ML LFSRs. 


A 



Consider the canonical representation of a rate (b + l)/b encoder [20] as shown in Fig. 2 when the 
switch is in position A. Let S k (D) be the state of the encoder at time k with coefficients , Sf, • ■ ■ , , 

where the output of the encoder at time k is 


x = s^_\ + Y y ^h t 


1 


(i) 


The state transition for input u k , • • • , u k at time k is given by 


S k {D) 


■ b 

Y, u i^(D) + DS k - l (D) 

.1=1 


mod h 0 (D) 


(2) 


From the all-zero state, we can enter the zero-input loop with nonzero input symbols tii, • • • ,tq, at state 
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( 3 ) 


b 

S l (D) = ' s ^2 / Uihi{D) mod ho(D) 

i=i 


From the^ same nonzero input symbol, we leave exactly at state S rn ~ 1 (D) back to the all-zero state, 
where S 2 > ~ 1 (D ) satisfies 


S l (D) = DS 2 ' n - l (D) mod h 0 (D) (4) 

i.e., S 2 ~ 1 (D) is the “predecessor” to state S 1 (D) in the zero-input loop. If the most significant bit of 
the predecessor state is zero, i.e., = 0, then the branch output for the transition from S 2 ' ~ l {D) 

to S l (D) is zero for a zero-input symbol. Now consider any weight-1 input symbol, i.e., Uj = 1 for j = i 
and Uj = 0 for j ^ z, j = 1, 2, ■ ■ • , b. The question is: What are the conditions on the coefficients h t (D) 
such that, if we enter with a weight-1 input symbol into the zero-input loop at state S X (D), the most 
significant bit of the “predecessor” state S r "~ l (D) is zero. Using Eqs. (3) and (4), we can establish that 

hio T hi m = 0 (5) 

Obviously, when we enter the zero-input loop from the all-zero state and when we leave this loop to go 
back to the all-zero state, we would like the parity output to be equal to 1. From Eqs. (1) and (5), we 
require 


hiQ — 1 
hi^rn ~ I 


> 


(6) 


With this condition, we can enter the zero-input loop with a weight-1 symbol at state S l (D) and then 
leave this loop from state 5 2 ~ 1 {D) back to the all-zero state, for the same weight-1 input. The parity 
w r eight of the codeword corresponding to weight-2 data sequences is then 2 m_1 -f 2, where the first term 
is the weight of the zero-input loop and the second term is due to the parity bit appearing when entering 
and leaving the loop. If b = 1, the proof is complete, and the condition to achieve the upper bound is 
given by Eq. (6). For b ~ 2, we may enter the zero-input loop with u = 10 at state S l (D) and leave the 
loop to the zero state with u = 01 at some state S j (D). If we can rhoose S j (D) such that the output 
weight of the zero-input loop from S l (D) to S J (D ) is exactly 2 m-1 /2 , then the output weight of the 
zero-input loop from S j+1 (D) to 5 2 ~ 1 (D) is exactly 2 m_1 /2, and the minimum weight of codewords 
corresponding to some weight-2 data sequences is 

2 m “ 1 


In general, for any 6, if we extend the procedure for 6 = 2, the minimum weight of the codewords 
corresponding to weight-2 data sequences is 


O 771— 1 

L— J + 2 (7) 

where [xj is the largest integer less than or equal to x. Clearly, this is the best achievable weight for the 
minimum-weight codeword corresponding to weight-2 data sequences. This upper bound can be achieved 


103 



if the maximum run length of l’s (m) in the zero-input loop does not exceed |_2 m 1 jb\. If m> [2 171 1 /6J , 
then the minimum weight of the codewords corresponding to weight-2 data sequences will be strictly less 
than [2 m - l /b\ +2. 

The run property of ML LFSRs [14] can help us in designing codes achieving this upper bound. 
Consider only runs of l’s with length l for 0 < / < m - 1; then there are 2 m ~ 2_/ runs of length /, no runs 
of length m - 1, and only one run of length m. ^ 

Corollary 1 . For any r — b/n recursive systematic convolutional code with b inputs, b systematic 
outputs, and n — b parity output bits using a primitive feedback generator, we have 


^2 - L 


(n - 6)2 


m — 1 


-J + 2(n - 6) 


(8) 


Proof. The total output weight of a zero-input loop due to parity bits is ( n - b)2 hl ~ 1 . In this zero- 
input loop, the largest minimum weight (due to parity bits) for entering and leaving the loop with any 
weight-1 input symbol is [(n — 6)2^ _1 ]/6. The output weight due to parity bits for entering and leaving 
the zero-input loop (both into and from the all-zero state) is 2(n - 6). & 

There is an advantage to using 6 > 1, since the bound in Eq. (8) for rate b/bn codes is larger than the 
bound for rate 1 fn codes. Examples of codes are found that meet the upper bound for b/bn codes. 

A. Best Rate b/b + 1 Constituent Codes 

We obtained the best rate 2/3 codes as shown in Table 1, where <^2 = ^2 + -• The minimum-weight 
codewords corresponding to weight-3 data sequences are denoted by ^3, is the minimum distance 

of the code, and k = m + 1 in all the tables. By “best” we mean only codes with a large (I 2 for a given 
m that result in a maximum effective free distance. We obtained the best rate 3/4 codes as shown in 
Table 2 and the best rate 4/5 codes as shown in Table 3. 


Table 1. Best rate 2/3 constituent codes. 


k 

Code generator 

d 2 

dz 

dniin 

3 

t- 

II 

0 

hi = 3 

h 2 = 5 

4 

3 

3 

4 

CO 

II 

0 

hi = 15 

h 2 = 17 

5 

4 

4 

5 

3- 

0 

1! 

to 

CO 

hi = 35 

h‘2 = 27 

8 

5 

5 


h 0 = 23 

hi = 35 

h 2 = 33 

8 

5 

5 

0 

h a = 45 

hi = 43 

h 2 = 61 

12 

6 

6 


Table 2. Best rate 3/4 constituent codes. 


k Code generator d 2 dz dmin 


ho - 7 

hi 

= 5 

h'2 = 3 

h 3 = 1 

3 

3 

3 

0 

II 

hi 

= 5 

h 2 = 3 

h 3 = 4 

3 

3 

3 

h 0 = 7 

hi 

Oft 5 

h 2 3 

hz = 2 

3 

3 

3 

h 0 = 13 

hi 

m 15 

h 2 = 17 

h 3 = 11 

4 

4 

4 

ho - 23 

hi 

= 35 

h 2 - 33 

hz = 25 

5 

4 

4 

h 0 = 23 

hi 

= 35 

r- 

CN 

II 

C4 

hz = 31 

5 

4 

4 

h 0 - 23 

hi 

= 35 

h 2 m 37 

hz - 21 

5 

4 

4 

h 0 = 23 

hi 

= 27 

h 2 = 37 

hz = 21 

5 

4 

4 
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Table 3. Best rate 4/5 constituent codes. 


k 


Code generator 


d-2 

d,3 

d’min 

4 

h 0 = 13 

hi = 15 

h 2 = 17 

h 3 = 11 

h 4 = 7 

4 

3 

3 


h 0 = 13 

hi = 15 

h 2 = 17 

h 3 = 11 

/14 = 5 

4 

3 

3 

5 

ho = 23 

hi = 35 

h 2 = 33 

h 3 = 37 

h 4 SB 31 

5 

4 

4 


h 0 = 23 

hi = 35 

h 2 = 27 

/i 3 = 37 

h 4 = 31 

5 

4 

4 


CO 

CN 

II 

O 

hi = 35 

h 2 = 21 

h 3 = 37 

/14 = 31 

5 

4 

4 



B. Trellis Termination for b/n Codes 

Trellis termination is performed (for b = 2, as an example) by setting the switch shown in Fig. 2 
in position B. The tap coefficients a^o, • * * , ai >m _ i for i = 1,2, ■ • ■ ,6 can be obtained by repeated use of 
Eq. (2) and by solving the resulting equations. The trellis can be terminated in state zero with at least 
m/b and at most m clock cycles. When Fig. 3 is extended to multiple input bits (b parallel feedback shift 
registers), a switch should be used for each input bit. 

C. Best Punctured Rate 1/2 Constituent Codes 

A rate 2/3 constituent code can be derived by puncturing the parity bit of a rate 1/2 recursive 
systematic convolutional code using, for example, a pattern P = [10]. A puncturing pattern P has zeros 
where parity bits are removed. 

Consider a rate 1/2 recursive systematic convolutional code (1, gi(D)/(g 0 (D)). For an input u(D), 
the parity output can be obtained as 


x(D) = 


niPMD) 

9o(D) 


( 9 ) 


We would like to puncture the output x(D) using, for example, the puncturing pattern P[10] (decimation 
by 2) and obtain the generator polynomials ho(D), hi(D), and / 12 (D) for the equivalent rate 2/3 code: 


1 0 


G = 


0 1 


MO)' 

/i 0 (D) 

/12(D) 

MD) . 
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We note that any polynomial f{D) = Yl a iD l , at e GF( 2), can be written as 


f(D) = h(D 2 ) + Df 2 (D 2 ) 


(10) 


where f\{D 2 ) corresponds to the even power terms of f(D), and Df 2 (D 2 ) corresponds to the odd power 
terms of f{D). Now, if we use this approach and apply it to the u(D), g\(D ), and go{D), then we can 
rewrite Eq. (9) as 


*1 (£> 2 ) + Dx 2 ( D 2 ) 


(«i {D 2 ) + Du 2 ( D 2 )) (gn {D 2 ) + D<? 12 (D 2 )) 
g<n(D*) + Dg 02 (D*) 


( 11 ) 


where x\ (D) and x 2 (D) correspond to the punctured output x(D) using puncturing patterns P[10] and 
P[01], respectively. If we multiply both sides of Eq. (11) by (goi(D 2 ) + P</o2(P 2 )) and equate the even 
and the odd power terms, we obtain two equations in two unknowns, namely X\(D) and x 2 (D). For 
example, solving for x\ (D), we obtain 


xi(D) = ui(D) 


hx(D) 

h 0 (D) 


+ u 2 (D) 


h 2 (P) 

ho(D) 


where ho(D) = go{D) and 


h\(D) = gn(D)goi(D) + Dg\ 2 (D)g{) 2 {D) 
h 2 (D) = Dg\ 2 (D)goi(D) 4- Dgn(D)go 2 (D) J 


(12) 


(13) 


From the second equation in Eq. (13), it is clear that /i 2 ,o = 0. A similar method can be used to show 
that for P [01] we get hi >rn = 0. These imply that the condition of Eq. (6) will be violated. Thus, wc have 
the following theorem. 

Theorem 2. If the parity puncturing pattern is P = [10] or P = [01], then it is impossible to achieve 
the upper bound on d 2 = d 1 ^ + 2 for rate 2/3 codes derived by puncturing rate 1/2 codes. 

The best rate 1 /2 constituent codes with puncturing pattern P = [10] that achieve the largest d 2 are 
given in Table 4. 


Table 4. Best rate 1/2 punctured 
constituent codes. 


k Code generator d 2 dmiii 


3 

go = 7 

91 

= 5 

4 

3 

3 

4 

go = 13 

9i 

= 15 

5 

4 

4 

5 

go = 23 

91 

- 37 

7 

4 

4 


g 0 = 23 

9i 

= 31 

7 

4 

4 


90 « 23 

91 

= 33 

6 

5 

5 


9o = 23 

9i 

= 35 

6 

4 

4 


9o = 23 

91 

m 27 

6 

4 

4 
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D. Best Rate 1/n Constituent Codes 


For rate 1/n codes, the upper bound in Eq. (7) for b = 1 reduces to 


d p 2 < (n - l)(2 m-1 +2) 


This upper bound was originally derived in [5], where the best rate 1/2 constituent codes meeting the 
bound were obtained. Here we present a simple proof based on our previous general result on rate b/n 
codes. Then we obtain the best rate 1 /3 and 1 /4 codes. 

Theorem 3. For rate 1/n recursive systematic convolutional codes with primitive feedback, we have 


d% < (n- l)(2 m_1 + 2) 


Proof. Consider a rate 1/n code, shown in Fig. 3. In this figure, go{D) is assumed to be a primitive 
polynomial. As discussed above, the output weight of the zero-input loop for parity bits is 2 771-1 inde- 
pendent of the choice of gi(D), i = 1, 2, • ■ • , n - 1, provided that g l {D) ^ 0 and that gt{D) j=- g 0 (D), by 
the shift-and-add and balance properties of ML LFSRs. If S(D) represents the state polynomial, then 
we can enter the zero-input loop only at state S l (D) = 1 and leave the loop to the all-zero state at state 
S 2 -1 (D) = D m ~ l . The ith parity output on the transition S 2 ’" ~ 1 (D) — ► S l (D) with a zero input bit is 


— 9i0 9i,m 


If 9 i 0 = 1 and g^ m = 1 for i = 1, • ■ • , n - 1, the output weight of the encoder for that transition is zero. 
The output weight due to the parity bits when entering and leaving the zero-input loop is (n — 1) for 
each case. In addition, the output weight of the zero-input loop will be (n - l)2 m_1 for (n - 1) parity 
bits. Thus, we established the upper bound on for rate 1/n codes. □ 


We obtained the best rate 1/3 and 1/4 codes without parity repetition, as shown in Tables 5 and 6, 
where d 2 = d% + 2 represents the minimum output weight given by weight-2 data sequences. The best 
rate 1/2 constituent codes are given by go and g\ in Table 5, as was also reported in [5]. 


Table 5. Best rate 1/3 constituent codes. 


k 

Code generator 

d 2 

^3 

d>min 

2 

<7o = 3 

9l = 2 

92 = 1 

4 

oo 

4 

3 

o 

II 

-4 

g\ = 5 

CO 

II 

cs 

cn 

8 

7 

7 

4 

go = 13 

91 = 17 

92 = 15 

14 

10 

10 

5 

90 = 23 

gi = 33 

CO 

II 

CM 

22 

12 

10 


<7o = 23 

91 = 25 

92 = 37 

22 

11 

11 
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Table 6. Best rate 1/4 constituent codes. 


k 


Code generator 


d-2 

d3 

dmin 

4 

9o = 13 

91 = 17 

92 = 15 

93 = 11 

20 

12 

12 

5 

go = 23 

9 i = 35 

92 = 27 

g 3 =37 

32 

16 

14 


go = 23 

9 \ = 33 

92=27 

93 = 37 

32 

16 

14 


go = 23 

9 i = 35 

92 = 33 

93 =37 

32 

16 

14 


go - 23 

9\ - 33 

92 =37 

93 = 25 

32 

15 

15 


E. Recursive Systematic Convolutional Codes With a Nonprimitive Feedback Polynomial 

So far, we assumed that the feedback polynomial for recursive systematic convolutional code is a 
primitive polynomial. We could ask whether it is possible to exceed the upper bound given in Theorem 1 
and Corollary 1 by using a nonprimitive polynomial. The answer is negative, thanks to a new theorem 
by Solomon W. Golomb (Appendix). 

Theorem 4. 1 For any rate 1/n linear recursive systematic convolutional code generated by a non- 
primitive feedback polynomial, the upper bound in Theorem 3 cannot be achieved, i.e., 


rig < (n- l)(2 m “ 1 +2) 


Proof. Using the results of Golomb (see the Appendix) for a nonprimitive feedback polynomial, there 
are more than two cycles (zero-input loops) in LFSR. The “zero cycle” has weight zero, and the weights 
of other cycles are nonzero. Thus, the weight of each cycle due to the results of the Appendix is strictly 
less than (n - l)2 m_1 . If we enter from the all-zero state with input weight-1 to one of the cycles of the 
shift register, then we have to leave the same cycle to the all-zero state with input weight- 1, as discussed 
in Theorem 1. Thus, d? < (n — l)(2 m “ 1 4- 2). Cl 

Theorem 5. For any rate b/b + 1 linear recursive systematic convolutional code generated by a 
nonprimitive feedback polynomial, the upper bound in Theorem 1 cannot be exceeded, i.e., 


om— 1 

dS < L-yJ + 2 

Proof. Again using the results of the Appendix, there is a “zero cycle” with weight zero and at least 
two cycles with nonzero weights, say q cycles with weights uq, w 2 , ■ * • , w q . The sum of the weights of all 
cycles is exactly 2 m_1 , i.e., w i — 2 m_1 . For a b/b + 1 code, we have b weight-1 symbols. Suppose that 
with b, of these weight-1 symbols we enter from the all-zero state to the ith cycle with weight then we 
have to leave the same cycle to the all-zero state with the same 6* symbols for i = 1, 2, ■ ■ ■ , <7, such that 
b L — b. Based on the discussion in the proof of Theorem 1, the largest achievable minimum output 
weight of codewords corresponding to weight- 2 sequences is mm(wi/b\, '^2/^2 > * * * , w q /b (I ) + 2. But it is 
easy to show that min('uq//>i, W2/b 2 , ■ ■ ■ , Wq/b q ) < Wi/Yl M = 2 m ~ l /b. A 


1 The proofs of Theorems 4 and 5 are based on a result by S. W. Golomb (see the Appendix), University of Southern 
California, Los Angeles, California, 1995. Theorem 4 and Corollary 2 were proved for more general cases when the 
code is generated by multiple LFSRs by R. J. McKliece, Communications Systems and Research Section, Jet Propulsion 
Laboratory, Pasadena, California, and California Institute of Technology, Pasadena, California, 1995, using a state-space 
approach. 
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Corollary 2. For any rate b/n linear recursive systematic convolutional code generated by a non- 
primitive feedback polynomial, the upper bound in Corollary 1 cannot be exceeded. 

Proof. The proof is similar to the Proof of Theorem 5, but now w i = ( n ~ b)2 rn ~ 1 . □ 


IV. Turbo Decoding for Multiple Codes 

In [9] we described a new turbo decoding scheme for q codes based on approximating the optimum 
bit decision rule. The scheme is based on solving a set of nonlinear equations given by (q = 3 is used to 
illustrate the concept) 


L 0 k 


= 2py 0 /c 


. = 1qct Eu:u,=i f(yi M II &k e ttj ^ i>,+ ^ 2;,+ ^ 3,) 

lfc ° S E u: u t =o ^(yil«) 

f , E u; u fc =i p (y2l u )n k e u ’( L '» +L 'j +L3 ^ 

Eu:ut =0 p ( y 2 |u) U^ k e^+^+L 3i ) 

£ = lo _ Eu^-1 p (y3l u ) ru 

^ " OS Eu:„ )k =0 P (y3|u)n j#fc ^(^+^+^) 

for k = 1, 2, * • • , N. In Eq. (14), Z^- represents extrinsic information and y*, i = 0, 1, 2, 3 ar e the received 
observation vectors corresponding to x*, i = 0, 1,2,3 (see Fig. 1), where p — y^rE^/TVo, if we assume 
the channel noise samples have unit variance per dimension. The final decision is then based on 


Lk = L ok + L ifc -F L2/C + L%k 


(15) 


which is passed through a hard limiter with a zero threshold. 

The above set of nonlinear equations is derived from the optimum bit decision rule, i.e., 


Lk = log 


Eu:u t =i p (yol u ) p (yil u ) p (y2l u ) p (y3|u) 

Eu:u k =o p (yo! u ) p (yi! ,, ) p (y2!u) p (y 3 !u) 


(16) 


using the following approximation: 


p (u|yi) ~ n 7 

k = 1 1 


N e u k Ljk 


+ e* 


(17) 


Note that, in general, P(u|y*) is not separable. The smaller the Kullback cross entropy [3,17] between 
right and left distributions in Eq. (17), the better is the approximation and, consequently, the closer is 
turbo decoding to the optimum bit decision. 
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We attempted to solve the nonlinear equations in Eq. (14) for L 1? L 2 , and L 3 by using the iterative 
procedure 


f(m+ 1 ) 
L \k 


= a[ m ^ log 


Eu:u t = l ^(yilu) n^ fc 


(18) 


for k = 1 , 2, • • - , N, iterating on m. Similar recursions hold for and L^\ The gain a ^ should 
be equal to one, but we noticed experimentally that better convergence can be obtained by optimizing 
this gain for each iteration, starting from a value less than 1 and increasing toward 1 with the iterations, 
as is often done in simulated annealing methods. We start the recursion with the initial condition 2 
= Lo- For the computation of Eq. (18), we use a modified MAP algorithm 3 with 
permuters (direct and inverse) where needed, as shown in Fig. 4. The MAP algorithm [ 1 ] always starts 
and ends at the all-zero state since we always terminate the trellis as described in [ 8 ]. We assumed 7T\ = I 
(identity); however, any 7 Ti can be used. The overall decoder is composed of block decoders connected 
as in Fig. 4, which can be implemented as a pipeline or by feedback. In [10] and [ 11 ], we proposed an 
alternative version of the above decoder that is more appropriate for use in turbo trellis-coded modulation, 
i.e., set Lo = 0 and consider yo as part of yi. If the systematic bits are distributed among encoders, we 
use the same distribution for y 0 among the MAP decoders. 



2 Note that the components of the Lj’s corresponding to the tail bits, i.e., for k — N + 1, ■ • • , N 4- Mi, are set to zero 
for all iterations. 

3 The modified MAP algorithm is described in S. Benedetto, D. Divsalar, G. Montorsi, and F. Pollara, “Soft-Output 
Decoding Algorithms in Iterative Decoding of Parallel Concatenated Convolutional Codes,” submitted to ICC ’ 96 . 
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At this point, further approximation for turbo decoding is possible if one term corresponding to a 
sequence u dominates other terms in the summation in the numerator and denominator of Eq. (18). 
Then the summations in Eq. (18) can be replaced by “maximum” operations with the same indices, i.e., 
replacing ^2 u:Uk=i with U '™V for i = 0 , 1 . A similar approximation can be used for and L^k in 
Eq. (14). This suboptimum decoder then corresponds to a turbo decoder that uses soft output Viterbi 
(SOVA)-type decoders rather than MAP decoders. Further approximations, i.e., replacing £ with max, 
can also be used in the MAP algorithm . 4 

A. Decoding Multiple Input Convolutional Codes 

If the rate b/n constituent code is not equivalent to a punctured rate 1 /n' code or if turbo trellis-coded 
modulation is used, we can first use the symbol MAP algorithm 5 to compute the log-likelihood ratio of 
a symbol u = U\,U 2 , ■ ■ ■ ,Ub given the observation y as 


A(u) = log 


Pi u|y) 
P( 0|y) 


where 0 corresponds to the all-zero symbol. Then we obtain the log-likelihood ratios of the jth bit within 
the symbol by 


L{uj) = log 


V e A(u ) 

Z^u:u, = l c 

y „e A < u > 


In this way, the turbo decoder operates on bits, and bit, rather than symbol, interleaving is used. 


V. Performance and Simulation Results 

The BER performance of these codes was evaluated by using transfer function bounds [4,6,12]. In [ 12 ], 
it was shown that transfer function bounds are very useful for SNR s above the cutoff rate threshold and 
that they cannot accurately predict performance in the region between cutoff rate and capacity. In this 
region, the performance was computed by simulation. 


Figure 5 shows the performance of turbo codes with m iterations and an interleaver size of N = 16, 384. 
The following codes are used as examples: 


( 1 ) Rate 1/2 Turbo Codes. 

Code A: Two 16-state, rate 2/3 constituent codes are used to construct a rate 1/2 turbo 
code as shown in Fig. 6 . The (worst-case) minimum codeword weights, corresponding 
to a weight-z input sequence for this code are d e j~ 14, ^ 3 = 7 , c? 4 = 8 , d 5 =5=d m ; n , and 
^ 6 — 6 . 


4 Ibid. 

5 Ibid. 
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Code B: A rate 1/2 turbo code also was constructed by using a differential encoder and a 
32-state, rate 1/2 code, as shown in Fig. 7. This is an example where the systematic bits 
for both encoders are not transmitted. The (worst-case) minimum codeword weights, 
corresponding to a weight-i input sequence for this code are d e f= 19, d 4 =6=d 7mn , de=9, 
d$= 8, and dio — 11- The output weights for odd i are large. 

(2) Rate 1/3 Turbo Code. 

Code C: Two 16-state, rate 1/2 constituent codes are used to construct a rate 1/3 turbo 
code as shown in Fig. 8. The (worst-case) minimum codeword weights, d*, corresponding 
to a weight-z input sequence for this code are d e f— 22, d 3 = 11, d 4 =12, d 5 = 9 = d mm > 
d 6 = 14, and d 7 =15. 

(3) Rate 1/4 Turbo Code. 

Code D: Two 16-state, rate 1/2 and rate 1/3 constituent codes are used to construct 
a rate 1/4 turbo code, as shown in Fig. 9, with d e f = 32, d 3 = 15 = d min , d 4 = 16, 
d 5 = 17, do — 16, and d 7 = 19. 

(4) Rate 1/15 Turbo Code. 

Code E: Two 16-state, rate 1/8 constituent codes are used to construct a rate 1/15 
turbo code, [l, g\/ 90 , 92 / 9o, 9s/ 90 , 9 a/ 90 , 95 / 9o,9o/ 9o,9i / 9o) and {g\/ 90 , 92 / 9o, 93 / 90 , 9a/ 
9o,9s/9o, 90 / 90 , 97 / 90 ), with go — (23) oc£a /, 9\ = (21 ) oc * a / , 92 — (25 ) oc tai, 9s — (27 ) oc tai, 
94 = (31 )octai, 95 - (33) octa/ , go = (35) octa *, and g 7 = (37) oc£a/ . The (worst-case) 
minimum codeword weights, d t , corresponding to a weight i input sequence for this code 
are d e /=142, d 3 =39=d m * n , d 4 =48, d 5 = 45, d 6 = 50, and d 7 =6 3. 

The simulation performance of other codes reported in this article is still in progress. 
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Fig. 6. Rate 1/2 turbo code constructed from two codes (/jq = 23, /7-j = 35, h 2 = 33). 



Fig. 7. Rate 1/2 turbo code constructed from a differential encoder and code 
(90 = 67, 9 1 = 73). 


VI. Turbo Trellis-Coded Modulation 

A pragmatic approach for turbo codes with multilevel modulation was proposed in [16]. Here we 
propose a different approach that outperforms the results in [16] when M-ary quadrature amplitude 
modulation (M-QAM) or M-ary phase shift keying (MPSK) modulation is used. A straightforward 
method for the use of turbo codes for multilevel modulation is first to select a rate b/(b + 1) constituent 
code, where the outputs are mapped to a 2 6+1 -level modulation based on Ungerboeck’s set partitioning 
method [21] (i.e., we can use Ungerboeck’s codes with feedback). If MPSK modulation is used, for every b 
bits at the input of the turbo encoder, we transmit two consecutive 2 b+1 phase-shift keying (PSK) signals, 
one per each encoder output. This results in a throughput of b/2 bits/s/Hz. If M-QAM modulation is 
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INPUT DATA 



Fig. 8. Rate 1/3 turbo code constructed from two identical codes 
(90 = 23, flf! = 33). 


INPUT DATA 



Fig. 9. Rate 1/4 turbo code constructed from two codes 
(SO = 23, Si = 33) and (go = 23, g<\ = 37, g^ - 25). 


used, we map the 6+1 outputs of the first component code to the 2 b+1 in-phase levels (I-channel) of a 
2" b+2 -QAM signal set and the 6+1 outputs of the second component code to the 2 b+l quadrature levels 
(Q-channel). The throughput of this system is 6 bits/s/Hz. 

First, we note that these methods require more levels of modulation than conventional trellis-coded 
modulation (TCM), which is not desirable in practice. Second, the input information sequences are used 
twice in the output modulation symbols, which also is not desirable. An obvious remedy is to puncture 
the output symbols of each trellis code and select the puncturing pattern such that the output symbols 
of the turbo code contain the input information only once. If the output symbols of the first encoder are 
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punctured, for example as 101010 • • •, the puncturing pattern of the second encoder must be nonuniform 
to guarantee that all information symbols are used, and it depends on the particular choice of interleaver. 
Now, for example, for 2 b+1 PSK, a throughput 6 can be achieved. This method has two drawbacks: It 
complicates the encoder and decoder, and the reliability of punctured symbols may not be fully estimated 
at the decoder. A better remedy, for rate 6/(6 + 1) (6 even) codes, is discussed in the next section. 


A. A New Method to Construct Turbo TCM 

For a q = 2 turbo code with rate 6/(6+ 1) constituent encoders, select the 6/2 systematic outputs and 
puncture the rest of the systematic outputs, but keep the parity bit of the 6/(6 + 1) code (note that the 
rate 6/(6 + 1) code may have been obtained already by puncturing a rate 1/2 code). Then do the same 
to the second constituent code, but select only those systematic bits that were punctured in the first 
encoder. This method requires at least two interleavers: The first interleaver permutes the bits selected 
by the first encoder and the second interleaver those punctured by the first encoder. For MPSK (or 
M-QAM), we can use 2 1+6/2 PSK symbols (or 2 1+6/2 QAM symbols) per encoder and achieve throughput 
6/2. For M-QAM, we can also use 2 1+6/2 levels in the I-channel and 2 1+6/2 levels in the Q-channel and 
achieve a throughput of 6 bits/s/Hz. These methods are equivalent to a multidimensional trellis-coded 
modulation scheme (in this case, two multilevel symbols per branch) that uses 2 6/2 x 2 1+6/2 symbols per 
branch, where the first symbol in the branch (which depends only on uncoded information) is punctured. 
Now, with these methods, the reliability of the punctured symbols can be fully estimated at the decoder. 
Obviously, the constituent codes for a given modulation should be redesigned based on the Euclidean 
distance. In this article, we give an example for 6 = 2 with 16-QAM modulation where, for simplicity, 
we can use the 2/3 codes in Table 1 with Gray code mapping. Note that this may result in suboptimum 
constituent codes for multilevel modulation. The turbo encoder with 16 QAM and two clock-cycle trellis 
termination is shown in Fig. 10. The BER performance of this code with the turbo decoding structure 
for two codes discussed in Section IV is given in Fig. 11. For permutations tt\ and 7 r 2 , we used S-random 
permutations [9] with S = 40 and S = 32, with a block size of 16,384 bits. For 8 PSK, we used two 
16-state, rate 4/5 codes given in Section V to achieve throughput 2. The parallel concatenated trellis 
codes with 8 PSK and two clock-cycle trellis termination is shown in Fig. 12. The BER performance of 
this code is given in Fig. 13. For 64 QAM, we used two 16-state, rate 4/5 codes given in Section V to 
achieve throughput 4. The parallel concatenated trellis codes with 64 QAM and two clock-cycle trellis 
termination is shown in Fig. 14. The BER performance of this code is given in Fig. 15. For permutations 
*1, 1 * 2 , 7T3, and 7T4 in Figs. 10, 12, and 14, we used random permutations, each with a block size of 4096 
bits. As was discussed above, there is no need to use four permutations; two permutations suffice, and 
they may even result in a better performance. Extension of the described method for construction of 
turbo TCM based on Euclidean distance is straightforward. 6 


VII. Conclusions 

In this article, we have shown that powerful turbo codes can be obtained if multiple constituent codes 
are used. We reviewed an iterative decoding method for multiple turbo codes by approximating the 
optimum bit decision rule. We obtained an upper bound on the effective free Euclidean distance of 6/n 
codes. We found the best rate 2/3, 3/4, 4/5, and 1/3 constituent codes that can be used in the design 
of multiple turbo codes. We proposed new schemes that can be used for power- and bandwidth-efficient 
turbo trellis-coded modulation. 


6 This is discussed in S. Benedetto, D. Divsalar, G. Montorsi, and F. Pollara, “Parallel Concatenated Trellis Coded Modu- 
lation,” submitted to ICC ’96. 
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Fig. 12. Parallel concatenated trellis-coded modulation, 8 P5K, 2 bits/s/Hz. 



Fig. 13. BER performance of parallel con- 
catenated trellis-coded modulation, 8 PSK, 
2 bits/s/Hx. 
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Appendix 

A Bound on the Weights of Shift Register Cycles 1 


I. Introduction 

A maximum-length linear shift register sequence — a pseudonoise (PN)-sequence or a maximal length 
(m)-sequence — of degree m has period p = 2 m - 1, with 2 m_1 ones and 2 m ~ 1 - 1 zeroes in each period. 
Thus, the weight of a PN cycle is 2 m_1 . From a linear shift register whose characteristic polynomial is 
reducible, or irreducible but not primitive, in addition to the “zero-cycle” of period 1, there are several 
other possible cycles, depending on the initial state of the register, and each of these cycles has a period 
less than 2 m - 1. 

The question is whether it is possible for any cycle, from any linear shift register of degree m, to have 
a weight greater than 2 m_1 . We shall show that the answer is “no” and that this result does not depend 
on the shift register being linear. 

II. The Main Result 

Let S be any feedback shift register of length m, linear or not. We need not even specify that the 
shift register produce “pure” cycles, without branches. We will use only the fact that each state of the 
shift register has a unique successor state. For any given initial state, we define the length L of the string 
starting from that state to be the number of states, counting from the initial state, prior to the second 
appearance of any state in the string. (In the case of branchless cycles, this is the length of the cycle with 
the given initial state.) 

The string itself is this succession of states of length L. The corresponding string sequence is the 
sequence of 0’s and l’s appearing in the right-most position of the register (or any other specific position 
of the register that has been agreed upon) as the string goes through its succession of L states. 

Theorem 1 . From a feedback shift register 5 of length m, the maximum number of l’s that can 
appear in any string sequence is 2 m ” 1 . 

Proof. There are 2 m possible states of the shift register S altogether. In any fixed position of the shift 
register, 2 rn_1 of these states have a 0 and 2 m_1 states have a 1. In a string of length L, all L of the states 
are distinct, and in any given position of the register, neither 0 nor 1 can occur more than 2 m “ 1 times. 
In particular, the weight of a string sequence from a register of length m cannot exceed 2 m_1 . □ 

Corollary 1 . No cycle from a feedback shift register of length m can have weight exceeding 2 m_1 . 


1 S. W. Golomb, personal communication, University of Southern California, Los Angeles, California, 1995. 


121 




N96- 16691 




6 3 



TDA Progress Report 42-123 


November 15, 1995 


The Trellis Complexity of Convolutional Codes 

R. J. McEliece 

Communications Systems and Research Section 
and 

California Institute of Technology 
Pasadena, California 

W. Lin 1 


It has long been known that convolutional codes have a natural , regular trellis 
structure that facilitates the implementation of Viterbi’s algorithm [30,10]. It has 
gradually become apparent that linear block codes also have a natural , though not 
in general a regular, “minimar trellis structure , which allows them to be decoded 
with a Viterbi-like algorithm [2,31,22,11,27,14,12,16,24,25,8,15]. In both cases, the 
complexity of the Viterbi decoding algorithm can be accurately estimated by the 
number of trellis edges per encoded bit. It would, therefore, appear that we are 
In a good position to make a fair comparison of the Viterbi decoding complex- 
ity of block and convolutional codes. Unfortunately, however, this comparison is 
somewhat muddled by the fact that some convolutional codes , the punctured con- 
volutional codes [4], are known to have trellis representations that are significantly 
less complex than the conventional trellis. In other words, the conventional trellis 
representation for a convolutional code may not be the minimal trellis representa- 
tion. Thus, ironically, at present we seem to know more about the minimal trellis 
representation for block than for convolutional codes. In this article, we provide a 
remedy, by developing a theory of minimal trellises for convolutional codes. (A sim- 
ilar theory has recently been given by Sidorenko and Zyablov [29].) This allows us 
to make a direct performance-complexity comparison for block and convolutional 
codes. A by-product of our work is an algorithm for choosing, from among all 
generator matrices for a given convolutional code, what we call a trellis-minimal 
generator matrix, from which the minimal trellis for the code can be directly con- 
structed. Another by-product is that, in the new theory, punctured convolutional 
codes no longer appear as a special class, but simply as high-rate convolutional 
codes whose trellis complexity is unexpectedly small. 


I. Introduction 

We begin with the standard definition of a convolutional code [9,26], always assuming that the under- 
lying field is F = GF( 2). An (n,k) convolutional code C is a A>dimensional subspace of F(D) n , where 
F(D) is the field of rational functions in the indeterminate D over the field F. The memory, or degree, 


1 Graduate student at the California Institute of Technology, Pasadena, California. 
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of C, is the smallest integer m such that C has an encoder requiring only m delay units. An (n, k) convo- 
lutional code with memory m is said to be a (n, k, m) convolutional code. The free distance of C is the 
minimum Hamming weight of any codeword in C. An (n, k, m) convolutional code with free distance d is 
said to be an (n, k, m, d) code. 

A minimal generator matrix G(D) for an (n, k, m) convolutional code C is a k x n matrix with polyno- 
mial entries, whose row space is C, such that the direct-form realization of an encoder for C based on G(D) 
uses exactly m delay elements [9,26]. From a minimal generator matrix G(D), or rather from a physical 
encoder built using G(D) as a blueprint, it is possible to construct a conventional trellis representation 
for C. This trellis is, in principle, infinite, but it has a very regular structure, consisting (after a short 
initial transient) of repeated copies of what we shall call the “trellis module” associated with G(D). The 
trellis module consists of 2 m initial states and 2 m final states, with each initial state being connected by 
a directed edge to exactly 2 fc final states. Thus, the trellis module has 2 fc+m edges. Each edge is labeled 
with an n-bit binary vector, namely, the output produced by the encoder in response to the given state 
transition. Thus, each edge has length (measured in edge labels) n, and so the total edge length of the 
conventional trellis module is n2 fc+m . Since each trellis module represents the encoder’s response to k 
input bits, we are led to define the conventional trellis complexity of the trellis module as 

” . 2 m+k edge labels per encoded bit (1) 

or edges per bit, for short. If the code C is decoded using Viterbi’s maximum-likelihood algorithm on the 
trellis [30,10], the work factor involved in updating the metrics and survivors at each trellis module is 
proportional to the edge length of the trellis module, so that the trellis complexity as defined in Eq. (1) 
is a measure of the effort per decoded bit required by Viterbi’s algorithm. (For a more detailed discussion 
of the complexity of Viterbi’s algorithm on a trellis, see [25, Section 2].) 

For example, consider the (3,2,2) convolutional code with minimal generator matrix given by 


Gi(D) 


(l + D 1 + D 1 \ 

\ D 0 1 + DJ 


(2) 


This code has the largest possible free distance, viz., d hee = 3, for any (3,2,2) code. A “direct-form” 
encoder based on the generator matrix G\{D) is shown in Fig. 1. If the input pair is (u\,U 2 ) and the 
state of the encoder is ( s,t ), then the output (xi,X 2 ,x 3 ) is given by 


Xi = Uj + s + t 


X 2 = Ui + s 


X 3 = Ui + u 2 +t 


(3) 


and the u next state” is just the input pair (t£i,u 2 ). The conventional trellis module for the code with 
minimal generator matrix G X (D) given in Eq. (2) is shown in Fig. 2. The three-bit edge label on the 
edge from (s,t) to (iti, u 2 ) is the triple (xi, 2 2 , x 3 ) given in Eq. (3). The total edge length is 48, so that 
the conventional trellis complexity corresponding to the matrix G X (D) is 48/2 = 24 edges per bit, as 
predicted by Eq. (1). 

But we can do substantially better than this, if we use the fact that this particular code is a punctured 
convolutional code. We now briefly review the theory of punctured convolutional codes to see how 
simplified trellises result. 
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Fig. 1. A direct-form encoder based on the generator 
matrix G-|(D) in Eq. (2). The input is (iq, U 2 ), the 
output is (x-j, X 2 , X 3 ), and the state of the encoder is 
(s, t). (The boxes labeled s and t are unit delay 
elements.) 



Fig. 2. The conventional trellis module 
for the code with minimal generator 
matrix Gi(D) given in Eq. (2). 


If we begin with a parent (N, l,?n) convolutional code, and block it to depth k, i.e., group the input 
bit stream into blocks of k bits each, the result is an ( Nk,k,m ) convolutional code. If we now delete, 
or puncture, all but n bits from each Nk - bit output block, the result is an ( n, fc,m ) convolutional code. 2 
This punctured code can be represented by a trellis whose trellis module is built from k copies of the 
trellis modules from the parent (N, l,m) code, each of which has only 2 m+1 edges, so that the total 
number of edge labels on the trellis module is n • 2 m+1 , which means that the trellis complexity of an 
(n,k,m) punctured code is 


~ • 2 m+1 edges per bit (4) 

k 

which is a factor of 2 k ~ 1 smaller than the complexity of the conventional trellis given in Eq. (1). For 
k — 1, this is no improvement, but for larger values of k, the decoding complexity reduction afforded 

2 In fact., the memory of the punctured code may be less than m, but for most interesting punctured codes, no memory 
reduction will take place. 


124 




by puncturing becomes increasingly significant. And while the class of punctured convolutional codes 
is considerably smaller than the class of unrestricted convolutional codes, nevertheless many punctured 
convolutional codes with good performance properties are known [4,13,3,7], and punctured convolutional 
codes, especially high-rate ones, are often preferred in practice. 

For example, consider the (2, 1,2,5) convolutional code defined by the minimal generator matrix 


G 2 (D) = (l + D + D 2 1 + D 2 ) (5) 

The conventional trellis module for this code is shown in Fig. 3. If we block this code into blocks of 
size k = 2, we obtain a (4,2,2) convolutional code, still with df ree — 5, for which the conventional trellis 
module is two copies of the trellis module shown in Fig. 3; see Fig. 4. 

Now we can do the puncturing. Take the (4, 2, 2) code, as represented by the trellis module in Fig. 4, 
and delete the second output bit on each of the edges in the second part of the module. The result is shown 
in Fig. 5. This structure can be thought of as the trellis module for a (3,2,2) code; the corresponding 
df ree turns out to be 3. According to Eq. (1), the conventional trellis complexity of a (3,2,2) code is 
3/2 • 2 4 = 24 edges per bit. But if we use instead the punctured trellis corresponding to the k — 2 
blocked version of the parent (2, 1,2) code, we find from Eq. (4), or Fig. 5, that the trellis complexity is 
instead only 3/2 • 2 3 = 12 edges per bit. In fact, it can be shown that this punctured (3, 2, 2) code is the 
same as the conventional code with generator matrix G\{D) given in Eq. (2). (Indeed, this example is 
taken almost verbatim from [4], where it was used to illustrate the way puncturing can reduce decoding 
complexity.) 





Fig. 3. The trellis module for 
the (2,i,2) code with genera- 
tor matrix G 2 (D) = (1 + D +D 2 
1 +D 2 ); total edge length is 
16, so the trellis complexity is 
16 edges per bit. 


Fig. 4. The trellis module for the (4,2,2) 
code obtained from the code of Fig. 3 
by blocking the inputs in blocks of size 
2; total edge length is 32, so the trellis 
complexity is 32/2 = 16 edges per bit, 
the same as for the original code. 


Fig. 5. The trellis module for the 
(3,2,2) punctured code obtained from 
the code of Fig. 4 by deleting every 
fourth bit; total edge length is 16 + 8 = 
24, so the trellis complexity is 24/2 = 
1 2 edges per bit. 


it seems mysterious that an ordinary-looking generator matrix like Gi{D) produces a code whose 
trellis complexity can be significantly reduced (if one knows that it is, in fact, a punctured code), whereas 
for an almost identical code, say one defined by the generator matrix 


(1 + D D 1 + D\ 

[Dll) 


no such reduction is apparently possible. In Section II, we will resolve this mystery by developing a 
simple algorithm for constructing the minimum possible trellis for any convolutional code. Our technique 
will always find a simplified trellis for a punctured code, with complexity at least as small as that given 
by Eq. (4), even if we are not told in advance that the code can be obtained by puncturing. But more 


125 


important, it will often result in considerable simplification of the trellis representation of a convolutional 
code that is not a punctured code. We will illustrate this with worked examples in Sections II and III 
and numerical tables in the Appendix. 


II. Construction of Minimal Trellises 

If G(D) is a minimal generator matrix for an (n, fc, m) convolutional code C, then we can write G(D) 
in the form 


G{D) = G 0 + G 1 D + --- + G l D l (6) 

where Gq," * ,Gl are k x n scalar matrices (i.e., matrices whose entries are from GF( 2)), and L is the 
maximum degree of any entry of G(D). If we concatenate the L + 1 matrices G o, * • • ,Gi, we obtain a 
k x (L+ 1 )n scalar matrix, which we denote by G: 


G = (G 0 G y G l ) 


( 7 ) 


It is well known [23, Chapter 9] that the matrix G and its shifts can be used to build a scalar generator 
matrix G sca i ar for the code C (for simplicity of notation, we illustrate the case L = 2): 



G 0 Gi 

g 2 




Go 

Gi 

g 2 


^scalar — 


G o 

Gi 

Go 

G 2 

G 1 G 2 


(8) 


The matrix in Eq. (8) is, except for the fact that it continues forever, the generator matrix for a binary 
block code (with a very regular structure), and so the techniques that have been developed for finding 
minimal trellises for block codes are useful for constructing trellis representations for convolutional codes. 
Here we apply the techniques developed in [25, Section 7], which show how to construct a trellis directly 
from any generator matrix for a given block code, and the minimal trellis if the generator is in minimal 
span form, to construct a trellis for C based on the infinite scalar generator matrix G sca i ar . 

The trellis module for the trellis associated with G sca i ar corresponds to the (L + l)k x n matrix module, 


( Gl \ 

G L - i 

Go ) 


( 9 ) 


which repeatedly appears as a vertical “slice” in G yca i ar . Using the techniques in [25, Section 7], it is easy 
to show that the number of edges in this trellis module is 


n 

edge count = ^^2 a ' (10) 

3 = i 
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where a.j is the number of active entries in the jth column of the matrix module G. (An element is called 
active if it belongs to the active span of one of the rows of G. We will elaborate on this below.) Our 
object, then, is to find a generator matrix for which the edge count in the corresponding trellis module is 
as small as possible. 

To clarify these ideas, we consider the (3,2, 1 ) code with (minimal) generator matrix 

g 3 (D) = (i 1+d l + D s j (11) 

According to Eq. ( 1 ), the conventional trellis complexity for this code is 12 edges per bit. However, we 
can do better. The scalar matrix G 3 corresponding to G 3 (D) is [cf. Eq. ( 7 )] 


O 3 


fl 0 1 0 0 0 \ 

V111011J 


(12) 


In Eq. (12), we have shown the active elements of each row, i.e., the entries from the first nonzero entry 
to the last nonzero entry, in boldface. The span length of (i.e., the number of active entries in) the first 
row_is, therefore, three; and the span length of the second row is six. The matrix module corresponding 
to G 3 is [cf. Eq. (9)] 

/0 0 0\ 

Oil 
10 1 
Vi 1 1 / 

Thus, ai = 3, <12 = 3, and 03 = 3, which by Eq. ( 10 ) means that the corresponding trellis module 
has 2 3 + 2 3 + 2 3 = 24 edges. Since each trellis module represents two encoded bits, the resulting trellis 
complexity is 24/2 = 12 edges per bit. Since we have already noted that the conventional trellis complexity 
for this code is also 12 edges per bit, the trellis corresponding to G 3 (£>) is not better than (in fact, it is 
isomorphic to) the conventional trellis. To do better, we need to find a generator matrix for the code for 
which V 2 a ' is less than 24. Using the results of [25, Section 6 ], it is possible to show that minimizing 
Et 2°' is equivalent to minimizing Et a *> Le -> the total span length of the corresponding G, and so we 
shall look for generator matrices for which the span of G is reduced. 

Note that if we add the first row of G 3 (D) to the second row, the resulting generator matrix, which is 
still minimal, is 



G »< D >-(° 1+ 0 O d) 


The scalar matrix G' 3 corresponding to G 3 (D) is [cf. Eq. ( 7 )] 


-ZlOlOOO 
3 ~ ^0 1 00 1 1 


(13) 


The span length of the first row of G' 3 (D) is three, and the span length of the second row is five, and so 

the total span length is eight, one less than that of G 3 (D). The matrix module corresponding to 51 is 
[cf. Eq. (9)] 
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G£ = 


/O 0 0\ 

Oil 
10 1 

\0 1 0 / 


Here a x =2 , a 2 = 3, and a 3 = 3, and so by Eq. (10) the corresponding trellis module has 2 2 + 2 3 + 2 3 = 20 
edges, so that the resulting trellis complexity is 20/2 = 10 edges per bit. The trellis module itself, 
constructed using the technique described in [25, Section 7] is shown in Fig. 0 . 



Fig. 6 . The trellis module for the ( 3 , 2 , 1 ) code with generator matrix 
G 3 (D). (Solid edges represent "0" code bits, and dashed edges 

represent “I" code bits. The labels on the vertices correspond to the 
information bits.) 


But we can do still better. If we multiply the first row of G' 3 (D) by D and add it to the second row, 
the resulting generator matrix, which is still minimal, is 


<%(D) 


1 0 1 \ 
D 1 + D OJ 


The scalar matrix G 3 corresponding to G 3 (D) is [cf. Eq. (7)] 


G’i = 


10 10 0 0 
0 10 110 


(14) 


The span length of G 3 (D) is seven, one less than that of G 3 (D). The matrix module corresponding to 
G'i is [cf. Eq. (9)] 

0 0 \ 

1 0 
0 1 
1 0 / 



Here ai — 2 , — 3, and a 3 — 2, and so by Eq. (10), the corresponding trellis module has 2 2 + 2 3 4- 2 2 = 16 

edges, so that, the resulting trellis complexity is 16/2 = 8 edges per bit. The trellis module itself, again 
constructed using the techniques described in [25, Section 7] is shown in Fig. 7. 
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Fig. 7. The trellis module for the (3,2,1) code with generator matrix 
G$(D). This is the minimal trellis module for this code. 


Furthermore, it is easy to see that there is no generator matrix for this code with span length less than 
seven, so that the trellis module shown in Fig. 7 yields the minimal trellis for the code. Alternatively, we 
examine the scalar generator matrix for the code corresponding to G'f [cf. Eq. (8)]: 


^scalar — 


‘I 0 1 

0 10 


0 0 0 
1 1 9 

1 o T 
o 1 o 


0 0 0 
1 I 0 
0 10 
0 1 0 


10 0 0 

l T o 


(15) 


In Eq. (15), we see that G sca i ar has the property that no column contains more than one underlined entry, 
the leftmost nonzero entry in its row (L), or more than one overlined entry, the rightmost nonzero entry 
in its row (R). Thus, G sca i ar has the LR property, and so, if it were a finite matrix, it would produce 
the minimal trellis for the code [25, Sec. 6]. To circumvent the problem that G sca iar is infinite, we can 
define the Mth truncation of the code C, denoted by as the ((M + L)n,Mk ) block code obtained 
by taking only the first Mk rows of Gscaiar i*e., the code with Mix (M -F L)n generator matrix 


■10 10 0 0 

0 1 0 1 I 0 

1 0 T 0 

nW = 0101 

scalar 


0 0 
I 0 


(16) 


10 10 0 0 

oio i To. 


Plainly, it G sca i ar has the Lk property, so does G^j , for all 


theory of trellises for block codes that the matrix G^ a } ar produces the minimal trellis for for all M, 
and so we can safely call the infinite trellis, built from trellis modules corresponding to G, the minimal 
trellis for the code. (Note that, in this example, the ratio of the conventional trellis complexity to the 
minimal trellis complexity is 12/8 = 3/2. If this code were punctured, then according to Eqs. (1) and 
(4), the ratio would be at least 2. Thus, we conclude that the code with generator matrix G 3 (D) as given 
in Eq. (11) is not a punctured code, which shows that the theory of minimal trellises for convolutional 
codes goes beyond merely “explaining” punctured codes.) 


The preceding argument, though it was presented in terms of a specific example, is entirely general. It 
shows that a basic generator matrix G(D) produces a minimal trellis if and only if G(D) has the property 
that the span length of the corresponding G cannot be reduced by an operation of the form 
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9i (D) ^ gi (D) + D e 9j (D) 


where g % {D) is the zth row of G(D) and £ is an integer in the range 0 < £ < L. We shall call a generator 
matrix with this property a trellis-minimal generator matrix for C. A trellis-minimal generator matrix 
must be minimal, but the converse need not be true, as the example of this section shows. Furthermore, 
it can be shown that the set of trellis-minimal generator matrices for a given code C coincides with the 
set of generator matrices for which the span length of the corresponding G is a minimum. In the next 
section, we will give two more examples of minimal trellises. 


III. Two More Examples 

Our first example is for the code whose generator matrix is given in Eq. (2). The corresponding 
decomposition [cf. Eq. (6)] is 


G,(D) 





D 


The scalar matrix G is, thus, 


G x = 



1 

0 


1 

1 


1 1 0\ 

1 0 l) 


and the matrix module G from Eq. (9) is then 


Gx 


( 1 1 °\ 

10 1 

111 

\0 0 1 / 


(17) 


Since there are three active entries in each column of G, it follows from Eq. (10) that the edge count for 
the trellis module is 2 3 -b 2 3 + 2 3 = 24, so that the trellis complexity for this trellis module is 24/2 = 12 
edges per bit, the same as given by Eq. (4) for the punctured trellis. To actually construct the trellis 
module, we can use the techniques of [25, Section 7], and the result is shown in Fig. 8. Finally, we note 
that the G sca iar corresponding to the matrix G\ of Eq. (17) is [cf. Eq. (8)] 


‘111110 
0 0 1 1 0 I 

1111 
0 0 11 
1 
0 


1 9 

0 1 

11110 
0 110 1 


which has the LR property, and so Gi(D) is trellis-minimal. (This code is the first code listed in Table 2 
in the Appendix.) 
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Fig. 8. The trellis module for the (3,2,2) code with generator matrix 
G-j(D). This module is isomorphic to the one in Fig. 5. 


As our second example, we consider a partial-unit-memory code, taken from [20,1], It is an (8,4,3) 
code with d free = 8 and with minimal generator matrix (as taken from [1]) 


G(D) 


/liiimix 
moiooo ’ 
10110100 
V 10011010/ 


oooooooo \ 
11011000 
10101100 
V 10010110/ 


(18) 


The conventional trellis complexity for this code is, by Eq. (1), 8/4 ■ 2 7 = 256 edges per bit. We can 
reduce this number to 120, as follows. First, we concatenate the two matrices in Eq. (18), obtaining the 
following 4 x 16 scalar matrix G : 


/I 11111110000000 0\ 
fll 1 010001 10110001 

I 1 0 1 1010010101 100 
\1 00110101001011 0 / 


Next, using the techniques developed in (25, Section 6j, we perform a series of elementary row operations 
on G, transforming it to the minimal span, or trellis oriented form, G'\ 


/I 11111110000000 0\ 
0001011111011000 ) 
0100101110101100 
\0 01011100011101 0 / 


(19) 


The matrix module G defined in Eq. (9) is, thus, 
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/O 0 0 0 0 0 0 0\ 

11011000 1 
10101100 
00111010 
11111111 
00010111 
01001011 
Vo 0 1 0 1 1 1 0/ 

and so by Eq. (10) the total edge length of the trellis module is 2 4 4- 2 5 + 2° 4- 2 7 4- 2 7 + 2 r> + 2 5 4- 2 4 = 480. 
Since each trellis module represents four encoded bits, it follows that the trellis complexity is 480/4 = 120 
edges per bit, compared to the conventional trellis complexity, cited above, of 250 edges per bit. 

The matrix G sca i ar corresponding to the matrix & in Eq. (19) is easily seen to have the LR property, 
and so the generator matrix [cf. Eq. (19)] 


/I 

1 

1 

1 

1 

1 

1 

1 \ 


(° 

0 

0 

0 

0 

0 

0 

°\ 

0 

0 

0 

1 

0 

1 

1 

1 

+ 

1 

1 

0 

1 

1 

0 

0 

0 

0 

1 

0 

0 

1 

0 

1 

1 

1 

0 

1 

0 

1 

1 

0 

0 

Vo 

0 

1 

0 

1 

1 

1 

0 / 


Vo 

0 

1 

1 

1 

0 

1 

0 / 


is trellis-minimal. However, the trellis complexity can be reduced still further, if we allow column per- 
mutations of the original generator matrix G(D) in Eq. (18). Indeed, by computer search, we have found 
that one minimal complexity column permutation for this particular code is the permutation (01243567), 
which results in the generator matrix [cf. Eq. (18)] 


G(D) = 


/ 11111111 \ 
11110000 1 
10101100 
V 10011010/ 


/ 00000000 \ 
11011000 ' 
10110100 
V 10001110/ 


( 20 ) 


Then, after putting the minimal generator matrix of Eq. (20) into trellis-minimal form, it becomes 


G(D) = 


/ 11111111 \ 
00001111 1 
01111111 
V 00111111/ 


/ 00000000 \ 
11111000 | 
11111100 
V 11111110/ 


D 


(21) 


The trellis complexity of the generator matrix in Eq. (21) turns out to be 104 edges per encoded bit. 
(This code is the seventh code listed in Table 6 in the Appendix.) The minimal trellis complexity of unit 
memory and partial unit memory convolutional codes has also been studied in [6] and [32]. 


IV. LTC Versus ACG 

In this section, we will attempt to compare the trellis complexity of a number of codes to their perfor- 
mance. To do this, we define the logarithmic trellis complexity (LTC) of a code, block or convolutional, 
as the base-2 logarithm of the minimal trellis complexity (edges per encoded bit) and the asymptotic 
coding gain (ACG) as the code’s rate times its minimum (or free) distance. An empirical study, based on 
existing tables of convolutional codes (e.g., the tables in [19,28,20,5,7]), reveals the interesting fact that 
LTC /ACG lies between 1.5 and 2.0 for most “good” convolutional codes. For example, for the (3, 2, 2, 3) 


132 



code discussed in Section III, the ratio is 1.79, and for the (8,4, 3, 8) code, it is 1.68. By comparison, for 
the “NASA standard” (2,1,6,10) convolutional code, for which, as for all (n, l,m) convolutional codes, 
the minimal trellis complexity is given by the formula of Eq. (1), the ratio is 1.60. In the Appendix, we 
list the (ACG,LTC) pairs for a large number of convolutional codes and a few block codes. In Fig. 9, 
we show a scatter plot of these pairs. It is interesting to note how close most of these pairs are to the 
line of slope 2. This experimental fact may be related to a recent theorem of Lafourcade and Vardy [18], 
which implies that for any sequence of block codes with a fixed rate R > 0 and fixed value of d/n > 0, as 
n — > oo, 


lim inf 

n— * oo 


LTC 

ACG 


> 2 


( 22 ) 


In any case, we have been able to show that for all codes, the ratio LTC /ACG must be strictly greater 
than 1. (This result is similar to Theorem 3 in [17].) 



Fig. 9. A scatter plot of the pairs (ACG, LTC) for the codes listed 
in the Appendix. 


V. Conclusion and Open Problems 

In this article, we have shown that every convolutional code has a unique minimal trellis representation, 
which is in many cases considerably simpler than the conventional trellis for the code. We have also 
presented a simple technique for actually constructing the minimal trellis for any convolutional code, 
and we have numerically computed the trellis complexity for many convolutional codes. In principle, the 
theory of minimal trellises for convolutional codes can be deduced from the general Forney-Trott theory 
[12], but we believe the observation that the Viterbi decoding complexity of many convolutional codes, 
including many nonpunctured codes, can thereby be reduced systematically is new, as are the details of 
the algorithms for producing the minimal trellises. 


We close with a list of research problems that suggest themselves. 
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(1) A given convolutional code will, in general, have many different minimal generator ma- 
trices [21], but as we saw in Section II, not all minimal generator matrices are trellis 
minimal. What can be said about the class of trellis minimal generator matrices? 

(2) A theoretical explanation of the experimental observation that most of the codes shown 
in Fig. 9 lie near the line of slope 2 would be welcome. 

(3) The design and implementation of Viterbi’s algorithm on conventional trellises is well 
understood. Since the techniques described here lead to greatly reduced trellis complex- 
ity, it would be worthwhile to make a careful study of how best to implement Viterbi’s 
algorithm on minimal trellises. 

(4) From our current viewpoint, punctured convolutional codes are just codes whose trellis 
module has fewer edges than would normally be expected. Indeed, it is easy to prove 
that the minimal trellis complexity of any punctured convolutional code is at least as 
small as the punctured trellis complexity given in Eq. (4). This is because in the scalar 
matrix G for a punctured code, certain entries are guaranteed to be zero. For example, 
for a (4,3,3) punctured code, the matrix G has the template structure 

( xxxxxxOO 
OQxxxxxO 
OOOxxxxx 

where the x’s can be arbitrary (actually, there are restrictions on the x’s that depend in 
detail on how the code is constructed), but the eight zero positions must be respected. 
Any (4, 3, 3) convolutional code with such a template structure will have trellis complexity 
at most 4/3 • 2 4 = 211/3. An obvious question is whether other low complexity templates 
support good convolutional codes. 

(5) In our computer-aided search for the “best” column permutation of the (8,4, 3,8) code, 
we found that each of the 8! = 40,326 possible column permutations had minimal trellis 
complexity of either 120 or 104. This strongly suggests an equivalence among permu- 
tations that, if understood theoretically, could make it much simpler to find the best 
column permutation. 

Finally, we remark that when the bulk of this article was written, we were not aware of the important 
earlier work of Sidorenko and Zyablov [29], which deals explicitly with the minimal trellis for a convo- 
lutional code, and we wish to acknowledge their priority. Their work, like ours, develops the theory of 
minimal trellises for convolutional codes from the corresponding theory for block codes. However, their 
trellis construction is based on the parity-check matrix of the code rather than the generator matrix, 
and their emphasis is quite different. One advantage of the Sidorenko-Zyablov approach is that it leads 
to the following upper bound on the number of nodes at depth i in the minimal trellis for a (n, k, rn) 
convolutional code [29, Theorem 1]: 


N t < 2 m+,nin ( /c ’ n_fc ) 


It is not easy to derive this bound using our methods. On the other hand, the present article contains a 
number of things not present in [29], among them being 
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(1) The observation that the minimal trellis for a punctured convolutional code is at least 
as simple as the punctured trellis. 

(2) The concept of a trellis- minimal generator matrix for a convolutional code, and an algo- 
rithm for computing one. 

(3) The ACG versus LTC comparison for block and convolutional codes. 
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Appendix 

Tables of LTC Versus ACG 


In this appendix, we list the ACG and the LTC for a large number of “good” convolutional codes and 
a few block codes. A scatter plot of these (ACG, LTC) pairs appears as Fig. 9 in Section IV. 


Table 1. Best (2,1 ,m) codes. a 


Code 

LTC 

ACG 

LTC- ACG 
ratio 

(2, 1,2,5) 

4 

2.5 

1.60 

(2, 1,3,6) 

5 

3 

1.67 

(2, 1,4,7) 

6 

3.5 

1.71 

(2, 1,5,8) 

7 

4 

1.75 

(2,1,6,10) 

8 

5 

1.60 

(2,1,8,12) 

10 

6 

1.67 

(2,1,10,14) 

12 

7 

1.71 

(2,1,11,15) 

13 

7.5 

1.73 

(2,1,12,16) 

14 

8 

1.75 

(2,1,14,18) 

16 

9 

1.78 

(2,1,15,19) 

17 

9.5 

1.79 

(2,1,16,20) 

18 

10 

1.80 

(2,1,18,22) 

20 

11 

1.82 

(2,1,21,24) 

23 

12 

1.92 

(2,1,23,26) 

25 

13 

1.92 

(2,1,25,27) 

27 

13.5 

2.00 

(2,1,27,28) 

29 

14 

2.07 

(2,1,30,30) 

32 

15 

2.13 

a From pp. 85-88 

in [7], 
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Table 2. Best (3,2, m) codes. 


Code 

LTC 

ACG 

LTC-ACG 

ratio 

(3, 2, 2, 3) 

3.58 

2.00 

1.79 

(3, 2, 3, 4) 

5.00 

2.67 

1.87 

(3, 2, 4, 5) 

6.00 

3.33 

1.80 

(3, 2, 5, 6) 

7.00 

4.00 

1.75 

(3, 2, 6, 7) 

8.00 

4.67 

1.71 

(3, 2, 7, 8) 

9.00 

5.33 

1.69 

(3, 2, 8, 8) 

10.00 

5.33 

1.88 

(3, 2,9, 9) 

11.00 

6.00 

1.83 

(3,2,10,10) 

12.00 

6.67 

1.80 


a From p 90 in [7]. 


Table 3. Best (4,3,m) codes. a 


Code 

LTC 

ACG 

LTC-ACG 

ratio 

(4, 3, 3, 4) 

5.00 

3.00 

1.67 

(4, 3,5, 5) 

7.00 

3.75 

1.87 

(4, 3, 6, 6) 

8.00 

4.50 

1.78 

(4, 3, 8, 7) 

10.00 

5.25 

1.90 

(4, 3, 9, 8) 

11.00 

6.00 

1.83 

a From p. 

90 in [7], 




Table 4. Best (3,1,m) codes. a 


Code 

LTC 

ACG 

LTC-ACG 

ratio 

(3, 1,2,8) 

4.58 

2.67 

1.72 

(3,1,3,10) 

5.58 

3.33 

1.68 

(3,1,4,12) 

6.58 

4.00 

1.64 

(3,1,5,13) 

7.58 

4.33 

1.75 

(3,1,6,15) 

8.58 

5.00 

1.72 

(3,1,7,16) 

9.58 

5.33 

1.80 

(3,1,8,18) 

10.58 

6.00 

1.76 

(3,1,9,20) 

11.58 

6.67 

1.74 

(3,1,10,22) 

12.58 

7.33 

1.72 

(3,1,11,24) 

13.58 

8.00 

1.70 

(3,1,12,24) 

14.58 

8.00 

1.82 

(3,1,13,26) 

15.58 

8.67 

1.80 


From p. 89 in [7j. 



Table 5. Best (4,1, m) codes/ 


Code 

LTC 

ACG 

LTC-ACG 

ratio 

(4,1,2,10) 

5.00 

2.50 

2.00 

(4,1,3,13) 

6.00 

3.25 

1.85 

(4,1,4,16) 

7.00 

4.00 

1.75 

(4,1,5,18) 

8.00 

4.50 

1.78 

(4,1,6,20) 

9.00 

5.00 

1.80 

(4,1,7,22) 

10.00 

5.50 

1.82 

(4,1,8,24 

11.00 

6.00 

1.83 

(4,1,9,27) 

12.00 

6.75 

1.78 

(4,1,10,29) 

13.00 

7.25 

1.79 

(4,1,11,32) 

14.00 

8.00 

1.75 

(4,1,12,33) 

15.00 

8.25 

1.82 

(4,1,13,36) 

16.00 

9.00 

1.78 

a From p. 89 

in [71. 




Table 6. Some block codes and partial unit memory 
convolutional codes. 

Code 

LTC 

ACG 

LTC-ACG 

ratio 

[8,4,4] 

Self-dual code 

3.46 

2.00 

1.73 

[24,12,8] 
Golay code 

8.22 

4.00 

2.06 

[32,16,8] 
BCH a code 

8.64 

4.00 

2.16 

[48,24,12] 
Self-dual code 

15.13 

6.00 

2.52 

[n, n - 1, 2] 
Parity-check code 

2.00 

2(n— 1) 
n 

n 

n— 1 

[n, 1, n] 

Repetition code 

1 + log 2 n 

1 

1 + log 2 Tl 

(8, 4, 3, 8) 
PUM b code 

6.70 

4.00 

1.68 

(24,12,7,12) 
PUM code 

15.58 

6.00 

2.60 

(24,12,10,16) 
PUM code 

18.58 

8.00 

2.32 

a Bose-Chaudhuri-Hocquenghem. 
b Partial unit memory. 
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This article describes measurements made at all three Deep Space Network 70-m 
S-band polarization diverse (SPD) systems to determine and eliminate the cause 
of the 1-K elevation in follow-up noise temperature in the listen-only mode of the 
SPD systems at DSS 43 and DSS 63. The system noise temperatures obtained after 
finding and correcting the cause of the elevated follow-up noise temperature are also 
reported . 


I. Introduction 

In response to the Galileo spacecraft’s X-band (8.45 GHz) antenna deployment failure, an emergency 
effort to optimize S-band (2.3 GHz) downlink performance was conducted. As part of this effort, termed 
the Galileo S-band Contingency Mission, the three 70-m DSN S-band polarization diverse (SPD) systems 
in the listen-only mode (see Fig. 1) have been carefully evaluated. Results of this initial evaluation were 
that both DSS 43 and DSS 63 at the Canberra and Madrid Deep Space Communications Complexes, 
respectively, exhibited elevated follow-up noise temperature contributions — defined as the contribution to 
system operating noise temperature of all components following the first low- noise amplifier (LNA) of 
1.25 K in comparison with the predicted values of 0.35 K. The system noise temperature predictions 
for these systems are shown in Tables 1 and 2. During the course of the evaluation, DSS-43 personnel 
determined the cause of this elevated follow-up noise temperature contribution in their 70-m SPD system 
to be due to a nonstandard configuration. This problem was corrected, and the antenna subsequently 
performed within predicted performance limits. 


The reason for the elevated follow-up noise temperature at DSS 63, also determined during the course 
of this work, was a higher than normal attenuation level in the system path behind the maser LNA. Once 
this attenuation was reduced by about 4 dB, the measured follow-up noise temperature at DSS 63 of 
0.4 K agreed very closely with the predicted value of 0.35 K. Also documented during the course of the 
investigation were high and dissimilar noise figures of the S-band Block IV receivers at DSS Tl, as well 
as differences in gain of the right-hand circular polarization (RCP) and left-hand circular polarization 
(LCP) channels of the very long baseline interferometry (VLBI) downconverter at DSS 63. 
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Fig. 1. The 70-m SPD diverse front-end listen-only mode for DSS 43 and 

DSS 63. 


II. Noise Measurement Technique 

In the course of this investigation, two types of noise measurements made at different points along the 
front-end component string were needed to assess where the noise temperature problem was located. The 
first measurement is the total system noise temperature, T op , while the antenna is pointed at zenith. The 
noise power level is measured using the 50-MHz precision attenuator assembly (PAS) or other suitable 
receiver. Switch 2 in Fig. 1 is switched so that the maser LNA is on the sky. Switch 2 is then switched 
so that the maser input is on the ambient load. The difference in power levels, Y, the Y-factor in dB, is 
then used to determine the T op : 


L op 


Tload + T ri 
10*71° 


(1) 


where Ti oa( i is the ambient load temperature, K, and T rcvr is the maser input noise temperature, K, 
including the follow-up noise temperature (FNT) contribution. This approach is extremely accurate for 
T rcv r « Tioad and requires independent knowledge of T rcvr . 

The second measurement required is the FNT. This is accomplished by measuring the difference in 
power level in dB, the Y-factor, by using the PAS or other suitable receiver, with the maser switched to 
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Table 1. Noise budget for DSS 14. a 


Component 

Part number 

AT, K 

L, dB 

T, K 

G, dB 

Tin , K 

Noise term, K 

Cosmic background 

- 

2.7 

— 

— 

— - 

— 

2.657 

Atmosphere 

— 

1.8 

— 

— 

— 


1.772 

Antenna spillover 

— 

1.44 

— 

— 

— 

— 

1.417 

Antenna scatter 

— 

2.3 

— 

— 

— 

— 

2.264 

Main reflector 


0.08 

— 

— 


— 

0.079 

Subreflector 

— 

0.07 

— 

— 

— 

— 

0.069 

Main reflector 

— 

0.1 

— 

— 

— 

— 

0.098 

gap leakage 

Feed horn 

9449420-1 

— 

0.003 

293 

— 

— 

0.199 

Waveguide round 

9457310-1 

— 

0.0015 

293 

-- 

— 

0.100 

Waveguide round, 

— 

— 

0.0018 

293 

— - 

— 

0.120 

15.141 in. 

+ rotary joint 1 in. 

Rotary joints (2) 

9457311-1 

— 

0.003 

293 

— 

— 

0.200 

Polarizers (1) 

9449405-1 

— 

0.0035 

293 

— 

— 

0.233 

Cosine taper 

9457389-1 

— 

0.002 

293 

— 

— 

0.133 

Orthomode, upper 

9457308-1 

— 

0.005 

293 

— 

— 

0.333 

Matched coupler, 

9457331-1 

— 

0.00345 

293 

— 

— 

0.230 

40 dB injected = 
0.029 K 

Elbow, H-plane 

9451160-2 

— 

0.0037 

293 

— 

— 

0.247 

S-band passband filter 

9430960 

— 

0.021 

293 

— 

— 

1.406 

3-position switch 

9443100-1 

— 

0.008 

293 

— 

— 

0.538 

Elbow, E-plane 

9451159-2 

— 

0.0037 

293 

— 

— 

0.249 

Straight, 13 in. 

9459426-3 

— 

0.003 

293 

— 

— 

0.202 

35-dB coupler (loss) 

SR8148D 

— 

0.0066 

293 

— 

— 

0.445 

35-dB coupler (injected) 

— 

— 

35 

293 

— 

— 

0.093 

Maser/CCR VSWR b 

— 

0.00 

— 

— 

— 

— 

0.000 

(not used) 

Maser/CCR package 

— 

— 

— 

— 

45 

2 

2.000 

LP filter 

— 

— 

0.1 

293 

— 

— 

0.000 

10-dB couplers (2) 

— 

— 

1 

293 

— 

— 

0.002 

Cable loss, 

— 

— 

1.7 

293 

— 

— 

0.006 

1/2 in. spiraline 

Maser select box 

— 

— 

0.8 

293 

— 

— 

0.004 

Cable 

— 

— 

0.5 

293 

— 


0.003 

Loss 

— 

— 

0 

0 

— 

— 

0.000 

Avantek amplifier 

AN-2200N 

— 

— 

— 

25 

870 

0.071 

20-dB coupler 

— 

— 

20 

293 

— 

— 

0.007 

Cable 

— 

— 

1 

293 

— 

— 

0.002 

Downconvcrter 7\ n 

— 

— 

— 

— 

— 

8881 

0.287 


Total antenna system noise temperature (referred to input of maser) 15.47 

Follow-up noise contribution Tf 0.382 


a SPD feedcone, low-noise path, 2295 MHz, 90-deg elevation angle, and clear weather. 
h Closed-cycle refrigerator (CCR) voltage standing-wave ratio (VSWR). 
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Table 2. Noise budget for DSS 43 and DSS 63. a 


Component 

Part number 

AT, K 

L, dB 

T, K 

G, dB 

Tin , K 

Noise term, K 

Cosmic background 

— 

2.7 

— 

— 

— 


2.656 

Atmosphere 

— 

1.9 

— 

— 

— 

— 

1.869 

Antenna spillover 

— 

1.44 

— 

— 

— 

— 

1.416 

Antenna scatter 

— 

2.3 

— 

— 

— 

— 

2.262 

Main reflector 

— 

0.08 

- 

— 

— 

— 

0.079 

Subreflector 

— 

0.07 

— 

— 

— 

— 

0.069 

Main reflector 

— 

0.1 

— . 

— 

— 

— 

0.098 

gap leakage 

Feedhorn 

9449420-1 

— 

0.003 

293 

— 

— 

0.199 

Waveguide round 

9457310-1 


0.0015 

293 

— 

— 

0.100 

Rotary joints (3) 

9457311-1 

— 

0.0045 

293 

— 

— 

0.299 

Polarizers (2) 

9449405-1 

— 

0.007 

293 

— 

— 

0.466 

Cosine taper 

9457389-1 

— 

0.002 

293 

— 

— 

0.133 

Orthomode, upper 

9457308-1 

— 

0.005 

293 

— 


0.333 

Matching section, 

9457331-1 

— 

0.003 

293 

— 

— 

0.200 

upper 

Elbow, H-plane 

9451160-2 

— 

0.0037 

293 

— 

— 

0.247 

S-band passband filter 

9430960 

— 

0.021 

293 

— 

— 

1.406 

3-position switch 

9443100-1 

— 

0.008 

293 

— 

— 

0.538 

Elbow, E-plane 

9451159-2 

— 

0.0037 

293 

— 

— 

0.249 

Straight, 13 in. 

9459426-3 

— 

0.003 

293 

— 

— 

0.202 

35-dB coupler (loss) 

SR8148D 

— 

0.0066 

293 

— 

— 

0.445 

35-dB coupler (injected) 

— 

— 

35 

293 

— 

— 

0.093 

Maser/CCR b package 

— 

— 

— 

— 

45 

2 

2.000 

LP filter 

— 

— 

0.1 

293 

— 

— 

0.000 

10-dB couplers (2) 

— 

— 

1 

293 

— 

— 

0.002 

Cable loss, 

— 

— 

1.7 

293 

— 

— 

0.006 

1/2 in. spiraline 

Maser select box 

— 

— 

0.8 

293 

— 

— 

0.004 

Cable 

— 

— 

0.5 

293 

— 

— 

0.003 

Avantek amplifier 

AN-2200N 

— 

— 

— 

25 

870 

0.071 

20-dB coupler 

— 

— 

20 

293 

— 

— 

0.007 

Cable 

- 

— 

1 

293 

— 

— 

0.002 

Downconverter Ti n 

— 

— 

— 

- 

— 

8881 

0.287 


Total antenna system noise temperature (referred to input of maser) 15.74 

Follow-up noise contribution Tj 0.382 


a SPD feedcone, low-noise path, 2295 MHz, 90-deg elevation angle, and clear weather. 
b Closed-cycle refrigerator. 
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the ambient load while switching the maser pump source on and off. The Y-factor, Y, is then used to 
determine the FNT: 


FNT 


Tload 


( 2 ) 


where Ti oa d is the ambient load temperature, K, and the difference in power level between the maser 
pump on and off is Y, dB. 


III. Preliminary Investigation and Baseline Data 

A noise budget was prepared for the DSS-14, DSS-43, and DSS-63 SPD systems in the listen-only 
mode. These noise budgets used our best estimates of microwave performance for each component in 
the system. Some measured data were available; other figures are theoretical. Measurements made at 
the stations were compared with these noise budget predictions. While the DSS-14 SPD system noise 
temperature agreed closely with its noise budget, those at DSS 43 and DSS 63 did not agree with predicted 
performance. Further FNT measurements isolated the problem at DSS 43 and DSS 63 to that part of 
the SPD system following the maser. This was determined after comparing the over-l-K FNT measured 
at both stations to the 0.4-K predicted noise. Since DSS 14 was the only station that closely agreed with 
predictions, and since it was the most readily available SPD system, it was carefully evaluated and used 
as a baseline against which to compare the other two stations. 

A T op measurement made at DSS 14 using the PAS yielded the data shown in row 1 of Table 3. 
A similar T op measurement was made at the immediate output of the maser using the JPL total power 
radiometer (TPR); this yielded the data in row 2. The T op measured at the input to the multiport coupler 
assembly (MPCA) is shown in row 3. An FNT measurement made using the PAS gave the data in row 4, 
while an FNT measurement made using the JPL TPR is shown in row 5. Next, the Block IV receivers’ 
noise performances were checked using a Hewlett Packard (HP) 8970B noise figure meter. The resulting 
noise figure (NF) and gain information obtained at 2295 MHz is displayed in row 6 for receiver 1 and 
row 7 for receiver 2. 

The high noise figures and the difference in noise figures, 21.4 dB for receiver 1 and 17.3 dB for 
receiver 2, of the Block IV receivers were noted. This poor performance results in the T op being more 
than 0.5-K above what good engineering practice should provide. A further explanation of the problem 
and a proposed solution appear in the recommendations section. 


Table 3. Measurement data for DSS 14. 


Data no., 
type 

y, dB 

Load, deg C 

Trcvr-i K 

r op , K 

FNT, K 

NF, dB 

Gain, « 

1, PAS 

12.88 

18 

2.5 

15.1 

— 

— 

— 

2, TPR 

13.05 

18 

2.5 

14.5 

— 

— 

— 

3, MPCA 

13.05 

18 

2.5 

14.5 

— 

— 

— 

4. PAS 

25 

18 

— 

— 

0.92 


— 

5, TPR 

35 

18 

— 

— 

0.04 

— 


6, Receiver 1 

— 

— 

— 

— 

— 

21.4 

33 

7, Receiver 2 

— 

— 

— 

— 


17.3 

33 
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IV. DSS-43 Measurements 


Measurements made at DSS 43 at the start of this investigation gave system noise temperature values 
of 17.2 K. Measurements made later on in the course of the investigation by station personnel at DSS 43 
using the 50-MHz PAS resulted in the data in rows 1 and 3 of Table 4. These data were reduced by 
station personnel and, therefore, the raw data are unavailable. Measurements made at DSS 43 using the 
JPL TPR yielded the data in rows 2 and 4. 

Station personnel explained the difference in measurements before and after evaluation activities as 
follows: At DSS 43, there never was a problem with the actual system noise temperature, and the precision 
power monitor method of measurement was reporting the correct result. The station chose to publish a 
determination of the system noise temperature based on a Y-factor detector result. This result was in 
error because of a nonstandard configuration of the Y-factor detector assembly. The station corrected 
this and confirmed that results from the three different methods 1 of measuring system noise temperature 
at the station agreed within 0.5 K. 

Table 4. Measurement data for DSS 43. 


Data no., 
type 

Y, dB 

Load, deg C 

Trcv ri K 

r op , K 

FNT, K 

NF, dB 

Gain, dB 

1, PAS 

— 

— 

— 

14.7 

--- 

— 

— 

2, TPR 

13.06 

18.8 

2.5 

14.7 

— 

— 

— 

3, PAS 

— 

— 

— 

— 

0.4 

N/A 

N/A 

4, TPR 

39.23 

19 

— 

— 

0.039 

N/A 

N/A 


V. DSS-63 Measurements 

The system operating noise temperature, T op , of the SPD system at DSS 63 was measured using three 
methods: Row 1 of Table 5 shows the measurement result using the 50-MHz PAS; row 2 shows the results 
using the JPL TPR at the output of the maser LNA; and row 3 shows the results using the JPL TPR at 
the input to the multiport coupler assembly. An FNT measurement made using the PAS gave the data in 
row 4, and an FNT measurement made using the JPL TPR is shown in row 5. The input noise figure of 
the Block IV receivers was measured using an HP 8970B noise figure meter, and the resulting noise figure 
and gain information obtained at 2295 MHz is displayed in row 6 for receiver 1 and row 7 for receiver 2. 
The VLB I downconverter input noise figure and gain for the R.CP and LCP channels were measured using 
an automated test setup developed at JPL. These data are shown in Fig. 2, which compares the noise 
figures of the two channels, and Fig. 3, which compares the gains. 

Included in earlier measurements were data taken of waveguide (WG) components on the ground using 
the S-band test horn, the TPR,, and a Block IV S-band maser that was brought from JPL. These tests 
were inconclusive due to the interaction of the WG components; that is, the WG components must be 
tested as a system and not independently in order to obtain accurate results. 

Plotting system noise temperature and ambient load temperature as a function of time (Fig. 4) revealed 
that the highest noise temperatures were occurring at the warmest load temperatures (i.e., the warmest 
time of the day). This was determined to be due to inadequate air conditioning in the cone, which was 
causing elevated physical temperatures of the WG components. The air conditioning was improved by 
station personnel to bring the SPD cone physical temperature down so that it more closely matched 
DSS 14 and DSS 43; this improved the stability of the system noise temperature over time. 

1 The third method uses the open-loop VLBI/radio science downconverter output and a quality spectrum analyzer for the 
Y-factor measurement. 
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Table 5. Measurement data for DSS 63. 


Data no., 
type 

y, dB 

Load, deg C 

Trcvr , K 

T op , K 

FNT, K 

NF, dB 

Gain, dB 

1, PAS 

12.54 

24 

4.5 

16.8 

— 


— 

2, TPR 

12.8 

27.5 

4.5 

16.0 

— 


— 

3, MPCA 

12.5 

27 

4.5 

17.1 

— 

- 

— 

4, PAS 

23.8 

25.5 

— 

— 

1.25 


— 

5, TPR 

39 

27 

— 

- 

0.04 

— 

— 

6, Receiver 1 

— 

— 

— 

— 

— 

15.7 

27 

7, Receiver 2 

— 

— 

— 

— 

— 

15.9 

30 



FREQUENCY, MHz 


Fig. 2. Comparison of the noise figures for the RCP and 
LCP channels of the VLBI downconverter at DSS 63. 



FREQUENCY, MHz 


Fig. 3. Comparison of the gain for the RCP and LCP 
channels of the VLBI downconverter at DSS 63. 
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Fig. 4. Comparison of T op and load temperature versus 
the time of day. 


VI. Conclusions 

The result of the evaluation was an overall improvement of the S-band polarization diverse systems 
at all stations except DSS 14. DSS 14 showed a measurable increase in system noise temperature per- 
formance; this is most likely due to an elevated noise temperature in the S-band Block IV receivers at 
the station. Pre-evaluation results showed a T op of 14.7 K, 17.2 K, and 17.6 K at DSS 14, DSS 43, and 
DSS 63, respectively. Postevaluation results showed a T op of 15.15 K, 14.7 K, and 15.5 K at DSS 14, 
DSS 43, and DSS 63, respectively. 

This task was successful in achieving its goals. The data taken and the equipment and procedures 
developed will assist in future investigations of station system noise temperatures. 


VII. Recommendations 

It is recommended that tolerances be established for the 70-m SPD system T op and FNT contributions 
at the stations. When an out-of-tolerance T op is measured, the FNT and linearity should be checked, 
and troubleshooting should proceed to identify the problem. This could and should include use of the 
HP 8970B noise figure meter with an HP 346-type noise source for the purpose of noise figure and gain 
measurements of station equipment behind the LNA. 

The Stelzried spreadsheet for checking Y-factor linearity should be implemented at all DSN stations. 
The stations should also have the capability of measuring the LNA noise temperature on the ground and 
at the output of the LNA in the cone, using the same system the JPL Microwave Electronics Group 
uses — a calibrated horn (for ground tests), absorber load, and the JPL total power radiometer. 

The attenuators and strip chart recorders should be replaced with a precision power meter or spectrum 
analyzer capable of measuring Y-factor power ratios at 50 MHz with an accuracy of ±0.01 dB. 

The T op of the 70-m SPD systems can be reduced an additional 0.3 K by reducing the follow-up noise 
temperature contribution by a factor of 10, from 0.4 to <0.04 K. This can be achieved by replacing the 
existing S-band follow-up amplifier with a state-of-the-art amplifier, installing this postamplifier in front 
of any losses between the LNA and the downconverter, and replacing the downconverters at the stations 
with state-of-the-art downconverters having noise figures of 5 dB or less. 
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The RCP channel of the VLBI downconverter at DSS 63 should be investigated. It has 14-dB less gain 
than the LCP channel. The fact that the RCP channel has a very similar noise figure when compared 
with the LCP channel indicates that this problem is in the output of the downconverter. 

Finally, the two S-band Block IV downconverters at DSS 14 exhibit elevated noise temperatures. This 
is most likely due to an elevated loss in either the preselector filter or the mixer and should be investigated 
and corrected. 
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